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Linear, time(shift)-invariant systems ene been exhaustively studied and 
their properties and behaviour are well known. These systems form the 
foundation of all engineering and scientific disciplines. However, they represent 
only an approximation of reality. This fact, of course, does not diminish their 
utility. Mathematical models employing the assumptions of linearity and shift- 
invariance often provide results sufficiently accurate to be of practical use. Fora 
large class of systems, however, these assumptions cannot be justified and so 
alternate models must be used. Mathematical models capable of representing 
nonlinearities and methods for their identification from system measurements are 
the major ‘topics explored in this thesis. Due to the particular treatment of 
nonlinear models chosen, multidimensional linear system modelling is also 
investigated. The assumption of shift-invariance will not be relaxed. 

It will be useful to formulate a geometric framework in which to solve the 
nonlinear modelling problem. The motivation for this is simple. The transition 
from physical problems to geometric ones allows many diverse phenomena to be 
handled with a common set of mathematical tools. This relieves the burden of 
having to invent new mathematics for each new situation. Instead. the well 
understood language of geometry is used to tackle many classes of problems. 
Kron |Ref. 1: p. 197] states very clearly the rationale that allows the real 
phvsical problem to be converted to an equivalent geometric one in the following 


passage. 


A set of n equations with n variables (with time as a parameter) may represent 
either the performance of a dynamical system with n degrees of freedom or the 
motion of a point along a curve located in an n-dimensional hypothetical space 
and expressed along some frame of reference. 


The basic approach taken in this dissertation. with respect to the modelling 
of nonlinear systems. is to represent them as a linear combination of nonlinear 
functions of the data. This allows linear algorithms to be applied in the solution 


of the modelling problem. In the process of solving the nonlinear problem several 
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new multidimensional lattice structures are developed. 

The inputs and outputs of the unknown system will be treated as random 
signals (not in general gaussian.) The problem of transforming random processes 
into geometric quantities has two solutions. These are introduced here in order 
to avoid confusion later and also in part to justify the tensor formulation that is 
used in the sequel. 

A random process KX may be defined as the assignment of a function 
{x(t.w), t€T}, to every outcome, w in a sample space 9. Of interest in this 
dissertation is the case when T is a discrete and finite index set. 

One way of geometrically visualizing this random process is to consider a 
Hilbert (or inner product) space, S, (in general infinite dimensional) of random 
variables. that is the vectors or elements of the space are random variables. 
Fixing t=ty x(tgw) 1s a random variable and so is a vector in S. The random 
process, X. a time series of random variables, is a series of vectors in S. or a curve 
in S [Ref. 2: p. 27]. The components of the vectors comprising X are indexed by 
the parameter w. The required inner product on this space is defined in terms of 
the statistical correlation, ie; E{x(t),w)x(t.w)}. This approach has proved highly 
successful in many applications [Ref. 3]. 

An alternate approach is to consider that the random process X. consists of 
vectors in a function space. In this interpretation the random process is an 
ensemble of time functions. {x(t,u),t<¢T}, indexed by w. Each of these time 
functions (generally referred to as realizations) is a vector in the function space. 
There will exist a large {in general infinite) number of such vectors corresponding 
to each possible outcome. «<<. The components of the vectors are indexed by 
the parameter t. There is no need to define a metric on this space. Any 
expectations that are required must be calculated over the ensemble of vectors. 

This second approach will be the one that is followed throughout this work. 
It will lead to many interesting and novel] interpretations of known algorithms 
and also will be used to derive several significant new results. 

While vectors are sufficient to provide a complete characterization of 


discrete-time. one-dimensional linear systems. general nonlinear systems with 
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memory require the use of higher order geometric objects to obtain convenient 
descriptions. It is shown in this dissertation that a particular class of these 
geometric objects that extend vector concepts in a natural way and provide an 
ability to deal with nonlinearities in an organized fashion are tensors. 

The contraversial Sapir-Whorf hypothesis from linguistics [Ref. 4] states 
that the constructs of a language define the boundaries of thought. 
Mathematics is a legitimate language. It is a well defined set of rules used to 
communicate ideas. If the mathematics that is employed in the solution of a 
problem is constrained, then it is conceivable that certain solutions may not be 
arrived at, or even that the problem may remain unsolved. In Electrical 
Engineering, vector calculus and linear algebra are the major mathematical tools. 
They are adequate to explain such diverse phenomena as the propagation of 
electromagnetic waves or the behaviour of one dimensional linear systems. More 
complex problems have also been solved using this theory by forcing them to fit, 
but the notation can become awkward. Tensor analysis is a convenient 
mathematical framework in which to deal with nonlinear and multidimensional 
signal processing problems. It provides an algebra for manipulating objects of 
higher dimension than two, which is all that can easily be handled using linear 
algebra. Importantly. tensor algebra furnishes a system of notation which is 
powerful, yet compact. 

The Electrical Engineer’s experience is usually limited to ordinary. 
Euclidean geometries. Physicists, around the turn of the century, began to 
realize that other more complex types of geometries were equally valid and 
important. In fact. Einstein showed that the world we live in is neither 
Euclidean nor is it simply three dimensional. 

The arguments outlined above provide the motivation for this study of the 
utility of tensor concepts in Electrical Engineering, specifically in the area of 
discrete signal processing. Although this work examines but a fraction of the 
possible applications in this field, it proves that tensors can lead to useful results 
and that they warrant further consideration, particularly in problems involving 


spaces of higher dimensions. 
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A. HISTORICAL BAGCKGHOeCND 

Tensor analysis has evolved in this century. originating with two Italian 
mathematicians in 1900: Ricci and Levi-Civita. Many of the early contributions 
are due to Einstein who required tensor concepts in the development of his 
general theory of relativity. Recent books on the subject include Golab, Synge 
and Schild, and Young [Ref. 5,6,7]. 

Kron [Ref. 1] in the early 1930’s made use of tensor concepts in Electrical 
Engineering. He appears to be first to do so. His work was mainly concerned 
with the analysis and design of electrical networks and rotating machinery. Since 
that time there have been few papers that deal with tensors in the context of 
Electrical Engineering. 

Volterra [Ref. 8] laid the foundation for nonlinear system analysis in the late 
nineteenth century. He studied functionals or functions of functions. He 
proposed a series of increasing order functionals as an approximation to any other 
functional. Frechet [Ref. 9: p. 517] later showed that this series was a complete 
representation and converged uniformly. This series has since become known as 
the Volterra series. We shall study this series in detail in Chapter 3. 

The first application of the Volterra series to nonlinear systems was done by 
Norbert Wiener in the 1930’s. Wiener’ also made several other significant 
contributions to nonlinear theory, such as the introduction of the Wiener G- 
functionals |Ref. 10]. They posses the property of orthogonality when the system 
input is white Gaussian noise. The two theories (Volterra and Wiener) form the 
basis of almost all significant work to date on nonlinear systems. 

One of the first practical methods of system identification was proposed by 
Lee and Schetzen [Ref. 11]. Their method takes advantage of the orthogonality of 
the Weiner G-functionals by emploving a cross-correlation technique to identity 
system parameters. 

The study of discrete-time nonlinear systems has gained importance with the 
advent of the digital computer. It appears that the idea of a discrete Volterra 
series first appeared in the mid 1960's (see for example [Ref. 12].) The use of 


tensor techniques in the study of nonlinear systems has received little attention. 
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Sandor and Williamson [Ref. 13] made some use of them in the study of 
continuous systems. More recently Parker and Thomas [Ref. 14] proposed the 
idea of using tensor methods for the analysis of nonlinear discrete-time systems. 
Their techniques for system identification involved the use of deterministic input 
signals to extract system parameters. 

The Volterra series is non-recursive and so a discrete form cannot represent 
an infinite memory system. This is equivalent to trying to represent an infinite 
memory linear system with a finite length impulse response. This can pose 
implementation difficulties for systems with long memories. One possible solution 
is the use of a recursive model. There has until very recently been little written 
about this because of the difficulty in analysing system stability. Parker and 
Perry [Ref. 15] have proposed a discrete nonlinear ARMA (auto-regressive 
moving-average) model. however. no stability implications were considered. Also 
Parker, Mayoral and Thomas [Ref. 16] proposed an Adaptive Kalman Identifier | 
or RLS (Recursive Least Square) type algorithm for the identification of non- 
linear ARMA systems. Zarzycki and Dewilde |Ref. 17] and Zarzycki [Ref. 
18,19,20,21] have proposed a nonlinear lattice structure. Again the stability of the 
resulting models is not discussed. Some nonlinear systems are inherently 
recursive (eg: the phase locked loop) so that this remains an important area for 
research. 

Recently, several books dealing exclusively with nonlinear systems theory 
have been published. The book by Schetzen [Ref. 9] concerns itself with 
continuous systems. It provides a very thorough but readable development of the 
classical concepts. Also of interest is a short appendix outlining the history of 
nonlinear systems theory. A book by Rugh [Ref. 22]. is an important contribution 


as it includes discussions of discrete theory. 


PeeOissERTATION OVERVIEW 
Because the typical reader of this dissertation will not have a background 
which includes tensor calculus it was felt that a chapter covering some 


fundamental concepts should be included. This material was considered to be of 
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central importance to the work that followed so it remains as a chapter rather 
than being relegated to an appendix. Readers familiar with tensor concepts may 
wish to skip most of Chapter 2. although. a cursory look is recommended to 
ensure that the notation is clearly understood. 

Chapter 3 begins with a review of the traditional Volterra theory of nonlinear 
systems. Both continuous and discrete-time systems are discussed. The 
discrete-time tensor equivalent of the discrete Volterra series is deduced. An 
alternate nonlinear tensor formulation is presented along with a discussion of its 
relationship to the Volterra series. This alternate formulation will be used in 
most of the work that follows. 

Next, methods for the identification of model parameters are examined. A 
tensor equivalent of the normal or Weiner-Hopf equations is formulated. Several 
recursive in time algorithms are included as examples of the application of 
traditional linear modelling methods to the nonlinear tensor formulation. 
Nontrivial numerical simulation results are included. 

The advantages of using alternate coordinate systems are then investigated. 
It is shown that by proper choice of coordinate systems and input signals the 
identification process can be significantly simplified. The nonlinear tensor 
formulation is extended to include recursive type models. It is shown that the 
Yule-Walker equations have a tensor counterpart which can be solved for the 
model parameters. Several of the new results of this chapter have already been 
published [Ref. 23]. 

In Chapter 4 a review of modern lattice theory is presented. Although the 
results themselves are not new the approach is novel. Tensor concepts are used 
to derive the lattice filters presented by considering orthogonalizing coordinate 
transformations. (Generalized forms of the Levinson and Schur algorithms are 
also presented and proven. These important algorithms are well known in linear 
matrix theory {Ref. 3] and their generalization in tensor form is a significant 
result. 

Chapter 5 breaks new ground by applying the lattice theory of Chapter 4 to 


the problem of modelling two-dimensional data fields. Simplifications due to an 
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assumption of shift-invariance are studied. Several different configurations are 
considered. Simulation results are included to prove the validity of the theory. 
Some implementational aspects of the algorithms are considered. In particular a 
systolic array is deduced for one of the two-dimensional] lattice algorithms 
presented. This result demonstrates that the new algortihms are amenable to 
implementation in dedicated VLSI hardware. 

In Chapter 6 a nonlinear lattice is formulated, again based upon the theory 
presented in Chapter 4 and 5. This is a new result. The lattice structure 
proposed differs from those of previous researchers in that it is recursive, not only 
in time order, but also in nonlinear order. For example, one can obtain the 
optimal cubic model from a knowledge of the optimal quadratic model. Once 
again nontrivial simulation results are included. 

Chapter 7 is a summary of the new results presented in this dissertation. It 
draws conclusions about these results and outlines some important unanswered 
questions as possible topics of future research. 

Two appendices are included. 

Appendix A contains an alternate proof of Theorem 4.4. This proof uses the 
Hilbert space formulation described in this introduction. It is included for two 
reasons: first. to illustrate this alternate formulation and second. to provide 
additional insight into this theorem which forms the foundation of Chapters 5 
and 6. 

Appendix B contains listings of the FORTRAN programs used in the 


simulations presented in this thesis. 
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Il MATHEMATICAL BACKGROUND 


This chapter presents a brief overview of the mathematical tools that are 
used in the dissertation. It begins with a discussion of linear, bilinear and 
multilinear functionals, and it is shown that these can be represented by tensors. 
Some customary conventions which simplify notation are introduced, and some 
useful tensor operations are presented and discussed. This is meant to be an 
introduction to the subject of tensor analysis. Only those concepts that will be 
used in the remainder of the dissertation are presented. Some proofs are included 
to act as examples. Many others are not presented here. since the interested 
reader can find them in the references [Ref. Bye), The discussion assumes a 


thorough knowledge of linear algebra. 


A. LINEAR FUNCTIONALS 

We say that. V, is a vector space over a field of scalars. F. if the operations 
of scalar multiplication and vector addition are defined such that the axioms of a 
vector space are satisfied (see for example [Ref. 24.25].) 

Consider a vector space. V, over a field of scalars, F. The elements of V are 
called vectors and will be denoted by use of boldface type, viz T. If we restrict 
otirselves to spaces of finite dimension. N, then we may write a vector T as an 
N-tuple of components and denote the vector space by V*. The components of a 
vector T will be denoted by T’. where A =1....N. Writing a vector as a set of 
components implies the existence of a basis. We will denote a basis for V* as the 
set of vectors A = {a,.....a,}. Thus an arbitrary vector T in V* can be written 
as a linear combination of these basis vectors. viz.. 


s 
(iis Sey (2.1) 


A=1} 


In order to maintain generality it is not necessary to commit to any specific 


vector space at this point. Likewise we allow the basis to remain arbitrary. 
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The following definitions and theorems are presented essentially without 
proof. 
1. Definition 2.1: Linear Functional 
If Vv‘ is a vector space over a field of scalars F, then, a linear 
transformation, H (the reason for the boldface will become apparent shortly, (see 
eqn (2.6))), of V‘ into F, is known as a linear functional or linear form on V~. 


We can indicate this transformation by 
H(T)=c where c€F and Te V* (222) 


2. Theorem 2.1 
The set of all linear functionals on V‘ forms a vector space of the same 
dimension as V*. This space is known as the dual vector space and is indicated 
ba 
3. Theorem 2.2 
If VN is a vector space over the field of scalars F, with basis A = 
{a,,..., ay}, then the set of linear functionals A = {b!-.... bY}, (defined so that the 
\-th functional. b’*. operating on an arbitrary element of V*. say T, yields the A- 
th component of T. namely T* ). form a basis for V*. The defining property can 


be expressed mathematically as 


DiT)=T° forallA €1.....N (2-3) 
N 

where T = aia, 
A= 


We call the functional b* the A-th coordinate function since when applied 
to a vector T it yields the A-th coordinate. namely T’. The set of these linear 
functionals, b>. A< {1..... N} comprising A is known as the dual basis of A. 

It is interesting to note that this choice of basis for the dual space leads 


to the property that 


] whendA= yp 
i a. oe vee 
b (a,,) ci ia {; whenA yp ( 


To show this we proceed as follows: from (2.1) it follows that 
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b'(T) = p Ms, (2.54) 


Sb'(ra,) (2.5b) 


p=1 
Since b’ is a linear functional (2.5b) can be written as, 
N 
oy) = YT*b* (a,) (2.5c) 
p=l1 
However, from (2.3) it is known that b*(T) = T*. Thus (2.5c) implies (2.4). 
The existence of a basis for the dual space implies that any vector (linear 
functional) H in V“ may be written uniquely as a linear combination of the 


elements of the dual basis. Therefore. any linear functional can be represented 


uniquely by an N-tuple of components. Thus 


From (2.6) one can write 


N 
H(T) = OH,b(T) (2. 


A=] 


to 

~] 

ev) 
~~” 


Using (2.3), (2.7a) becomes 
N 
IsBb) = Oy 15.10” (2.7b) 
A=] 
Alternately. in matrix form. if 
eee 0s, oe! (238) 
T (2.9) 
then (2.6) can be written as: 
H(T) = HT. (2.10) 


We notice that the two vectors. H and T. are defined relative to different 


basis and that they belong to different vector spaces. The vector H. defined 
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according to the dual basis A is called a covariant vector. The vector T. defined 
according to the regular basis. A. is known as a contravariant vector. As a 
convention. whenever a subscript is used to index the components of a vector it 
is understood that the vector is being expressed according to the dual basis, A, in 
the dual vector space V*, and is a covariant vector. Similarly, when a 
superscripted index is used the components are assumed to represent a 
contravariant vector in the vector space V™ according to the regular basis A. 

Equation (2.7b) represents what we normally think of as a vector inner 
product. We usually do not think of the two vectors as coming from different 
vector spaces. The reason for this is that in the familiar rectangular cartesian 
system of coordinates, the regular and dual basis are identical and so there is no 
need to differentiate between covariant and contravariant vectors. In other 
systems of coordinates the distinction must be made in order that the 
relationships have meaning. To perform a vector inner product. one vector must 
be covariant and the other must be contravariant. For example consider the 
vector T as illustrated in Figure 2.1. It has components 1 37 with respect to 
basis {e,,e.} and components ‘-1 2)" with respect to basis {e,-.e,-}. 

A measure of the length of vector T in the rectangular coordinate system. 


{e;.e,}, can be computed from the expression 








3 1/2 
ar = | TT (2.11a) 
A=1 
9 *) 12 
7 ]12~ 2°| (2.11b) 
co: (2.11c) 


The answer is correct because in the rectangular Cartesian coordinate system 
there is no need to distinguish between covariant and contravariant vectors. le: 
T*= T,. However. an expression similar to (2.11c) in the oblique coordinate 


system, {e,..e,:}. does not yield a measure of the length of the vector. 
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Figure 2.1: Vector T expressed according to the two bases {e,.e,} and {e,-.e.-}. 


i) 
ty 


N 1/2 
ya = \(—1)2+2| (2.12a) 
A=1 


= 51? (2.12b) 


The answer is incorrect since both vectors in expression (2.12) are contravariant. 
To obtain a correct answer, one of the vectors would have to be made covariant. 
This involves introduction of a metric tensor. The interested reader can find a 
discussion of this concept in |Ref. 6]. In the case of rectangular coordinate 
systems this metric tensor is the identity matrix. Our intent here was only to 
indicate that in any system other than rectangular Cartesian, strict attention 


must be paid to the character of the vectors. 


B. TENSORS 

In the previous section the concept of covariant or contravariant vectors has 
been established. In this section a more formal definition of these quantities is 
presented. 


‘then a set of values of these 


Suppose we have N variables x’. x*.....x 
variables is called a point. The variables themselves are called coordinates (or 
components.) The totality of all points. as each of the variables (coordinates) 
x*, A = 1,.... N. vary over their entire specified range. constitutes an N-dimensional 
space, denoted by V*. 

1. Definition 2.2: Contravariant Vector 

A contravariant vector T. is defined on the basis of the transformation of 
its components upon transition from one coordinate system to another. For 


coordinate system (A) the components of T are an N-tuple of nuinbers designated 


as 
fe (A = J... N) 


Upon transition to another coordinate system {A‘). if the components of T 


transforni according to the rule 


23 


<1 (2.13) 


where x* and x* define the coordinates of a point in the old (A) and new ()’) 
coordinate systems, then T is said to be a contravariant vector. 
2. Definition 2.3: Covariant Vector 
A covariant vector U, is defined on the basis of the transformation of its 
components upon transition from one coordinate system to another. For 
coordinate system (A) the components of U are an N-tuple of numbers designated 


as 


Upon transition to another coordinate system (\’), if the components of U 


transform according to the rule 





xem 
Uy (2.14) 


> 
U,. = 
mie oe: 


then 1s Saidsto be a covariant vecion 


Onc 


In equation (2.13) the quantity = 
xX 


represents the partial derivative of 





the new (primed) coordinates with respect to the old coordinates. Similarly. in 


ax? 


el represents the partial derivative of the old coordinates with 
x 


equation (2.14) 





respect to the new. pritned. coordinates. In general these quantities can be 
arranged into a two-dimensional matrix of numbers. 
Ge looeug tole, 2. Ih 


Consider the following parametric description of a curve: 


Bes (2 Tea 
Se Gli (2. Tei 
aN iia) ( 2a 


We consider the three quantities x’, x*.x*, to be the coordinates of a point or 


equivalently the components of a vector in a three dimensional space. We leave 


the basis arbitrary. Indeed the equations we write will be true regardless of the 
choice of basis. This is the inherent advantage of tensor analysis since it allows 
expressions to be written which are invariant with respect to the coordinate 
system. | 

We can now define new, primed, coordinates as functions of the old 


coordinates. For example, if we arbitrarily choose the following coordinate 





transformation 
x =) e4(x',x*,x°) = \/ is (2.16a) 
a ceiKx. xX) = =x? (2.16b) 
moe = Cex xa) = —(2x*°—x') \ (2.16c) 
then, the quantities “— can be written in matrix form as 


tQ 
poe 
~] 





Ox?’ J+ 
= 0 —_— O 
-\/= , ~ / > 
TT 8 


According to equation 2.13 any contravariant vector. say T. with components T’. 





can be expressed in the new (primed) coordinate system as 





or in terms of components 


T= fer _ oT? - 0Tt - fer (2.182) 
vis a 

T? = oT! 4 fir + 0-T? = Jar (2.18b) 
Tv vfs 
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T =. /+t + OT? + 2m» / at fe ur — T’) (2, lege 
7 7 T 


4. Definition 2.4: Contravariant Tensor of Order 2 
A set of N? numbers T*“, where A and » = 1,...,.N are said to be the 
components of a second order contravariant tensor if, upon transition to another 


coordinate system, they transform according to the rule 








N N De 
eS S 7 om dx" (2.19) 


A=l1 p=l Ox Ox# 


5. Definition 2.5: Covariant Tensor of Order 2 
A set of N? numbers U,,. where A and » = 1,....N are said to be the 
components of a second order covariant tensor if, upon transition to another 


coordinate system. they transform according to the rule 


\ 








2.20) 


Similarly, tensors of higher order can also be defined. In the general 
case. it will no longer be possible to use different letters to denote indices. In this 
case indices with sub-indices will be utilized, namely Ay, A», .... An. 


o- 


6. Definition 2.6: Contravariant Tensor of Order p 
yee 
A set of NP numbers T"! >, where 4, = 1, ..., N fori = 1,...,p. are said to be 
the components of a p-th order contravariant tensor if. upon transition to 


another coordinate system. they transform according to the rule 





: ; N N ay a 
d d d 41 Ox Ox 
] p ioe 1 Pp tae 

i 7 (2.21) 

Ay } “A Ba ] ‘e] x ox Pp 

7. Definition 2.7: Covariant Tensor of Order p 
A set of NP numbers U. 4. WwhG@re 217g Sige eee ree sald to be 
‘D 


the components of a p-th order covariant tensor if. upon transition to another 


coordinate system. they transform according to the rule 








~~ (232) 


8. Definition 2.8: Mixed Tensor of Order p 


idteugae 
A set of NP? numbers S' “4, To Where. = 1. ..., Nforr= 1,...,.p. are Said 
Pp 


q+] 
to be the components of a p-th order mixed tensor. with q contravariant and (p- 
q) covariant indices if, upon transition to another coordinate system, they 


transform according to the rule 











fe ° 
1 a 
= A ql 2"? 
we ae d i 
= B ie As ox, pic Mae)» < ig aon? 
= S wae os ee (2.23) 
\,=1 A j=l ox OX EO co ° TO ” 


We have already seen two examples of mixed tensors. The first is the 


Kronecker delta 


1 whendA = yp 
a : 
a f whend + py (2 24) 


To see that the Kronecker delta is in fact a mixed tensor we must test to see if it 
transforms according to the rule given in equation (2.23). We must prove that 


the following relation is true 























N oN ie 
ax. Ox" 
il ore 25 
We begin with the right hand side of (2.25) 
NN c N 7 aad 
Or OK Ox xe 
—- -§* = 6° UY: 
Py py Ox? Ox# Py ax Py “axt | a) 
No eco ) 
ox Ox 
= 2G 
4 ox Oxe | ‘i 
ax? 
aoe (2206 |} 
ngs 3 
en (22ed) 


Therefore, we conclude that relation (2.25) holds and so the Kronecker delta is in 


fact a second order iensor of mixed character. 
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The other mixed character tensor that we have already worked with is 
the one appearing in the formulae of transformation (definitions 2.2 through 2.8). 
that is the the partial derivative of the new coordinates with respect to the old 
(and the old with respect to the new). We will not prove that this is in fact a 
tensor of the type stated although the proof is straight forward. It is instructive 


x anal Ox 
Ox? ax? 


are inverses of each 








to note, however, that the two quantities 


other. 
As a final note, vectors are tensors of order one. Also scalars are 
considered to be tensors of order zero. They are sometimes called invariants 


since their representation is independent of the coordinate systems used. 


C. BNOTATIONAL CONN ENTIONS 

There are two widely accepted conventions that simplify notation and 
unquestionably save much writing. The first is known as the summation 
convention. Historically, it was first used by Einstein. He noticed that in almost 
all cases there is really no need to explicitly write summation symbols. 
Summation can be implied whenever an index is repeated in an expression. once 
as a superscript and once as a subscript. The repeated index is allowed to take 
on all permissible values and the resulting terms are summed together. This type 
of index is often referred to as a dummy index. But what are permissible values 
for the index? This question leads to the second convention. Normally. the range 
of the index will not be explicitly stated. By convention it is understood that all 
greek subscripts and superscripts appearing in an expression will take on all 
values from 1 to N. where N is the dimension of the vector space in which we are 
working. In later chapters we will find it more convenient to allow indices to run 
froin O to N. The dimensionality of the vector space will thus be N+1. An 
additional convention which we shall find useful 1s to reserve latin indices to 
indicate that we are dealing with a particular component of a quantity. In most 
books this is indicated by surrounding the particular index with parenthesis. 
However. we will reserve parenthesis to indicate exponentiation. The conventions 


adopted here will be used throughout the sequel. In exceptional cases. where 
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some deviation from them is required. we will explicitly state the meaning of the 


notation. As an example of the use of these conventions consider the expression 
N 
eae for A = 1, ...,N D2 
p=) 


It can be written more succinctly as: 
Ya = yy” ; (2.28) 


Another convention which has already been used is now formally introduced. 
Every tensor quantity will be given a distinct base letter. Upon a change of 
coordinates, the base letter will be maintained in order to indicate that the 
quantity has not been modified, only the representation has changed. The 
coordinate system used is indicated by the sub- or superscript used to index the 
components of the quantity. We therefore, will refer to different coordinate 
systems simply by the index letter used to indicate the components. For 
example, a vector T has components T* in the (A) set of coordinates, while it has 
components T*’ in the (A’}) coordinate system. We note that using this 
representation, scalars appear identical in all coordinate systems. which is 
desirable. 

1. Example 2.2 

The following example serves not only as an illustration of the 
conventions presented in this section, but also as a concrete (and presumably a 
somewhat familiar) illustration of the two types of vector. Consider an invariant 


function of the coordinates. f = f(x’, x’. x°). The differential of this function is given 


by 


Of oO ‘one 
ene ee 8 
Ox Ox" Ox 


dx? (2.29] 


We can consider dx* to be the components of a contravariant vector representing 
an infinitesimal displacement expressed according to some basis A = {a,. ao. az}. 


Peeecan write this as: 
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N 
dx = )) dx’a, (2.30) 


A=1 


We have already shown that components of a contravariant vector transform 


according to 


se ch! (2.31) 


Ox? 





The gradient is also a vector whose components are given by: 


Of 


ax? 





Vf, = (2922) 


Upon a change of coordinates the values of the new components, f,-, can be 


deduced by application of the chain rule; 


A 
ae ppg galt (2.33) 


dx* Ox? 





It is clear that the components of the gradient transform according to the rule 
given in equation (2.14). The gradient must, therefore. be considered to be a 
covariant vector. 
2, Example 2.3 

Although it has already been stated (section 2.1) that linear functionals 
can be considered to be covariant vectors, this fact has not been proven. In this 
example we will show that any linear functional. say H. wit hy GOrpOnents H, thas 
transforms an arbitrary contravariant vector. say T. (according to equation 
(2.7b)) to yield an invariant. satisfies the definition of a covariant vector 
(equation (2.14)). Equation (2.7b). which defines a vector inner product is 


repeated here for convenience. 
Jal Ie) ee ed (2.3-4a} 


The quantities T* are the components of an arbitrary comtravariant vector. 


Because of the assumed invariance we mnay write 
Hl =e (2.34b) 


From the definition of a contravariant vector, (equation (2.13)) 
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ax? 














fe = 1" 5 [ 22546) 
Equation (2.34b) becomes 
ae 
ee = Hy Te (2.34d) 
Ox 
Rearranging yields, 
ax? >| 
ia — Hi, =) Sa (2.34e) 
Ox 
Equation (2.34e) implies 
Ox? 
H, = H,- 30 (2.34f) 


since the contravariant vector, T, was arbitrary. A simple change of variables in 


this last expression yields 


which is identical to the equation defining a covariant vector. (2.14). We are 


thus justified in calling linear functionals covariant vectors. 


i BiILINEAR AND MULTILINEAR FORMS 
We next consider higher order functionals. In particular we will start with 
the bilinear form or bilinear functional. 
1. Definition 2.9: Bilinear Functional 
A mapping f of a pair of vectors. say T<¢ U* and Se V™ into a field of 
scalars. is a bilinear functional or bilinear forin if f(T.S) is a linear function of 
T and S taken independently. We will only consider cases when N = M and 
= Vv. 
Once again choosing an arbitrary basis A = {a,. .. a,} we may express 


two arbitrary contravariant vectors as linear combinations of these basis vectors. 
ie Ta,, and S = S“a, 


The bilinear functional f can then be written 


ol 


f(S.1) = t(S> oe an (2.35a) 
= 5 fl ayeear) (2.35b} 
= 5° fiaan) (2.35) 


The bilinear form is thus completely determined by the N’ quantities f{a,,a,). We 
will write these components of f as f,,. Using this shorthand, equation (2.35c) can 


be written as 

[Sena oa (2.35¢) 
In matrix notation the bilinear form takes on the familiar appearance 

f(S.T) =sseF Tt (2.36) 


where F = !f,,] 
If the two vectors S and T are equal then the bilinear form reduces to the 


well known quadratic form 

a emt re a Be (2. 37am) 
or in Matrix notation 

{lt tre (2. 30iBy 


We will be interested in the behaviour of the components. f,,. of the 
bilinear form. f. upon transition from one coordinate system to another. We 
establish their tensor character in the following example. 

2. Example 2.4 

In Example 2.3 we showed that a linear functional satisfied the definition 
of a covariant vector. Here we will show that a bilinear functional satisfies the 
definition of a covariant tensor of second order. It is necessary for the discussion 
that follows in later chapters to establish the tensor character of bilinear 
functionals. Since the bilinear functional vields an invariant {scalar}. we may 


Byte 


f 39° Tn oe (2.38a) 


o2 


Ox? 


ee T 





a ae pe (2.38b) 
Ox" Ox? Xd 
= uy pel fi2yco 1 (2.38c} 
Rearranging yields 
Ox# Ox?” 
oe or buen )S*T =aG | (2.38d) 


Since the contravariant vectors S and T are arbitrary the quantity inside the 
parenthesis must vanish. This yields the relation 


Ox4# ax? 


Ox dx? 


(2.38e} 


iT? i pr 


This last equation is identical to the definition of a second order covariant tensor 
(equation (2.20).) We have thus proven that bilinear functionals are covariant 
tensors of order 2. 

In general we can have m-linear functionals which map m contravariant 
vectors. T(1),.... T(m), into a scalar and are linear functions of each of the m input 
vectors taken separately. They can be written as 


man. Tim))= T 1) --- T mf, . ., (2.39) 


1 m 


Using identical arguments to the ones presented for linear and_ bilinear 


functionals. we can show that multilinear functionals are also covariant tensors. 


fle NSOR OPERATIONS 

There are a few tensor operations that will be of considerable importance in 
later discussions. Although some have already been used we will formally define 
them here before proceeding. Only those operations that will be used in the 
sequel are presented. Others are possible and are discussed at length in the 


references (see for example [Refs. 1.5,6.7].) 
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1. Definition 2.10: Tensor Outer Product 
Given the components of two tensors 


T ly “Sh | i. and gAn+ Ae (2.40) 


the N®*9) numbers R' ae given by 


Poth a Anja A 
fen = ca. ? - (2.41) 


Him4+1 


are components of a tensor of order p + q. The operation implied in the above is 
known as the tensor outer product or simply the tensor product. 
2. Example 2.5 


As an example of the tensor product operation consider two vectors 


defined as: 

ae Send Se S| (2.42) 
The tensor product is given by (equation (2.41}) 

Rig eel Se (2.43) 


In this case the components of the tensor product can be arranged as a matrix of 


N* numbers. For simplicity consider the case when N=3. In matrix notation 


Re Ts: (2.44a) 
Or 
ROE eee (2.44b) 
T'S! T's? Tis? 
- Ip2%s! T2s?2 T2032 (2.44) 
T?s! T*S- T?S? 


3. | Definition 21 Conia om 

The operation of setting two indices. one lower and the other upper. 
equal and summing the result is known as contraction. The result 1s a tensor 
which has the character indicated by the remaining indices. The contraction of a 


tensor of order p + 1 over two indices results in a tensor of order p- 1. 


o4 


4. Example 2.6 


We may contract the tensor 


over the indices 4 and yz, resulting in 


H’, = H', + H*, + H’; (2.46a) 


II 


Trace|H?, | (2.46b) 


a scalar which represents the trace of the matrix. 
We can consider a higher order example. Assuming a three dimensional 
vector space. consider contracting a tensor U*,” over the indices \ and ». The 


result will be a vector whose components, U7, are given by 


mee eeu) + Ue, (2.47a) 
Wee! = U2,’ + UF,” (2.47b) 
Pe], = U*,” — v*,? [217¢) 


5. Definition 2.12: Tensor Inner Product 


Suppose that after we form the outer product of two arbitrary tensors T 


cone 3 ie er 
fees with components T' °°, ., ands"? . we set the indices 4. 


uy jin et H, 


and yw, equal. This implies a contraction operation. The outer product is 


written (according to equation (2.41)) as 


oo? a aie See, (2.48a) 


uy Hy uy u 


‘oa u y# 


f u Hy He yr fay i 


jo q 
It can be shown that the object. R. given in the last equation is a tensor of the 


character shown. If the original tensors. T and S. were of order m+n and p+q- 


(m+n) respectively. then the resulting tensor. R. will be of order (p+q-2). The 
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total operation. consisting of contracting the result of a tensor product over one 
(or more) pairs of indices. each member of the pair having come from a different 
tensor. is known as a tensor inner product. or simply inner product. This 
operation is also sometimes referred to as transvection (see for example [Ref. 
5].) 

This operation has been used previously in the discussion of linear, 
bilinear and multilinear functionals, (see equations (2.7b), (2.35d) and (2.39) 
respectively.) It will be particulariy valuable in future discussions. 

6. Example 2.7 

Some insight into the inner product operation may be gained by 

explicitly performing the two steps described above for the simple case of the 


linear functional (equation (2.7b).) 
ae (2.49) 


We can first form the outer product by replacing one of the 4 indices appearing 


on the right hand side with a different letter. say ». We may then write: 


Sg Wee aed ee 
HT) = eee (250) 
Hel 2 Hike 


The result is now contracted by equating A and » and performing the implied 


summation. 
y= int ae (251a) 
= ul 1 55 He a Het (2 ouniey 


It is often useful to perform both these steps (tensor outer product and 
contraction) when calculating an inner product. particularly when dealing with 
higher order tensors. It may otherwise be difficult to keep track of how terms are 


combined or even which terms will be present in the final expression. 


hie NONLINEAR SYSTEMS THEORY 


In this chapter several models of nonlinear systems are described. The 
discussion begins with an overview of the classical continuous Volterra series. 
Several interesting aspects of the series are examined. Next, a discrete-time 
version of the Volterra series is introduced and is used to develop an equivalent 
tensor form. Finally, a new discrete-time, nonlinear tensor model is presented 


and its relationship to the Volterra series is discussed. 


PeOCONTINUOUS NONLINEAR SYSTEMS MODELS 
Traditionally. non-linear systems have been modelled using the continuous 
Volterra series expansion |Ref. 26: p. 1559]. In its most general form this can be 


written as: 


y(t) = ho(t) 


=a i h,(t.7 ix(F ridin 


~ ff holtr ir e)x(r i)x(r 2)d7 1dr» 


= 00 oC 


<7 et 


where x(t} is the svstern input and y(t} is the system output. The parameters 
ho(t). hy(t.7,). ho(t.7,.72). °° are known as the Volterra Kernels. As we will only 
be concerned with time-invariant systems. the kernels will only be functions of 
the time difference. (t-7). and not the actual time. The kernels then becorne: 


fit —7 j). he{t-7,,t-7.). ---. Simple changes of variables then allow the Series. 


(3.1), to be written in the form 


O7 


aD f hte 1)x(t—7)dry 


+ ff holrara)x(t—7 y)x(t—r 2)d7 dr, 


— 00= 00 


res (3.2) 


We note several things about the expansion. First, an infinite number of 
terms are required to represent the most general case. Second, the kernels 
ho, hy(7 1), ha(r,,72), °°: correspond to the constant, linear, quadratic, ... terms of 
the expansion, respectively. The familiar linear system appears as a special case 
of this more general expansion (ie; the case when all kernels except h,(7,) are 


zero.) The third term: 


f hol T 1,7 9)x(t—7 ,)x(t—75)dr ,dr. (3.3) 


oF 10.6 
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a bilinear term. that is it is linear in each variable x(t—7,) and x(t-—7.) taken 
independently (ie: assume the other variable is a constant.) In general. the 
(i¢+1)-st term is i-linear. It is linear in each of the i variables 
x(t—r,). x(t-—7 9), .... x(t-7,) taken one at a time. Lastly. notice that the Volterra 
expansion is not orthogonal in the sense that the identification of the n-th order 
kernel depends on the values of all the other kernels. They cannot be identified 
independently. 

The non-linear model represented by equation (3.1) can be visualized as 
shown in Figure 3.1. As can be seen it corresponds to a parallel connection of 
subsystems of increasing non-linearity. Each of the subsystems is homogeneous 
(except for the zeroth order subsystem) in the sense that increasing (multiplying) 
the input by a factor k results in the output of the p-th subsystem being 
increased by a factor k?. This can easily be understood by examining the 


expression for the p-th term. 
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i : oi. i T »)x(t—7 4) eet adr ye > dr, (3.4) 


I “70a fara, Nt -7 1) °° > kx(t—7.Jdr, °*° dr, (3.5a) 
= a aa f hyo. taxi 7) 7° Xitoe dr, tee din [3255)) 


The presence of a constant term in the Volterra expansion should not be 


unexpected. Consider, for example, a system whose response is given by: 
y(t)=x(t}+h (3G) 


It is easily shown that this system does not obey the principle of superposition 
and so cannot be considered linear. This necessitates the inclusion of a constant 
term in the Volterra expansion in order to handle the general case. The usual 
procedure adopted in linear analysis. if a constant term appears, is to define a 
new output function which is the actual output less the constant term. This new 
output function is then identified in the usual fashion. The constant term must 
be separately identified. If we admit that the system is non-linear then this will 
no longer be necessary. 

There are many other aspects of continuous-time nonlinear system modelling 
that have not been discussed. In the next section discrete-time nonlinear systems 
are introduced. Many of the comments that will be made there are equally 
applicable to continuous-time nonlinear systems. However. we make them in the 


context of discrete-time, since that will make them more applicable to the sequel. 


meee I~CRETE-TIMENONDINEAR SYSTEM MODELS 


A discrete version of the time-invariant Volterra expansion can be written as: 


y(k} = ho 


ae > hy(Ay}x(k—A) 


A,=0 
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Figure 3.1: Nonlinear Volterra Series Model 
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er WAIL. )xtk —A, x(k - A) 


ies. 


—_ (3.7) 


where we have assumed that the system is causal. 

For a large class of systems the summations can be truncated to N+1 terms. 
Certainly truncation is required for computer implementation. This truncation 
implies that only finite memory nonlinear systems can be represented. Equation 
(3.7) is an extension of the linear Moving Average (MA) type model. In fact the 
linear model is a special case of equation (3.7). This expansion is non-recursive in 
nature, that is. it expands the present output only in terms of the present and 
past input. Past outputs are not used. 

The representation of the kernels in both the continuous and discrete forms 
of the Volterra expansion is not unique. There are. however, several special forms 
which are important. Consider a second order’ kernel for which 
ho(A,,A2) = he(Az.A,). This kernel is symmetric with respect to the two parameters 4, 
and \.. It turns out that the kernel can always be symmetrized with no loss of 
generality. For the p-th order kernel the procedure for obtaining the symmetric 


kernel from an asymmetric one is given by |Ref. 22]: 
CS , 
sym(Ay- «+» Ag) = Be Shale s+: Ag(p)) (3-8) 


where the summation is over all n! possible permutations of the p A’s. Although 
the svmmetric form may contain more terms than an asymmetric form it is of 
importance because it is unique |Ref. 9: p. 43]. There can be many equivalent 
asymmetric forms of the kernel which all lead to the same symmetric kernel 
(through equation (3.8})}. The symmetric kernel thus provides a standard forin 
Which can be used as a reference. 

iter forms of interest are also possible. The symmetry of the synimetric 
kernel imnplies redundancy. This redundancy can be eliminated by use of a 


triangular kernel. Consider a kernel! defined so that h,,,(A;. .... A,) = 0 Whenever 


4] 


\,>A, for s < t. The domain of a second order triangular kernel is illustrated in 
Figure 3.2a. For comparison. the equivalent symmetric kernel’s domain is 
illustrated in Figure 3.2b. 

Using the triangular kernel the output of the p-th order nonlinear subsystem 


is given by 


Ap As N 
y (k) oe yy » ‘ah s hei (Ay, os Ap)x(k— Ay) —- x(k AS) (Sane) 
A,=0 A,=0 Ag 


Notice that the limits of the summations reflect the triangular domain. This 
implies that fewer terms are included in the summations resulting in 
computational savings. The triangular kernel defined above, and used in 
equation (3.10). is not unique. Other triangular kernels can be formed by 
choosing alternate triangular domains. For example, the domain illustrated in 
Figure 3.3 could equivalently have been used. This choice corresponds to setting 
h n(AiseeA,) =@iwhenever Ay > wien = 7. 

The output of nonrecursive models of the type presented in equation (3.7) is 
stable if the input is bounded and if the series is truncated to a finite number of 
summations. Consider an input x(n) which is bounded to be less than some 
constant M. If the series is truncated to p+1 terms then. in the worst case the 
output will be 


N es N 
yin) <h-+ MS) hy(Ay)) + M? 3) DS E he(AyAz) + 
NiO A= 0 


A,=0 


ae) > es y h,(A4. ieee A.) (3.8) 


So. as long as the summations are bounded (which will generally be the case). the 
output will always remain finite. This guaranteed stability makes MA type 
models very attractive. As mentioned earlier. their shortcoming is their inability 
to accurately model infinite memory systems without using a large number of 


TeTIUS. 


42 


N 


a. Second Order Triangular Kernel Domain 





b. Second Order Syminetric Kernel Domain 


Figure 3.2: Triangular and Symmetric Volterra Kerne! Domains 


1. Tensor Formulation 
In order to adopt a tensor formulation of the problem we notice that 
equation (3.7) can be considered to consist of a series of increasing order 
functionals. As has been shown, these functionals can be expressed as tensors. 


We can therefore rewrite equation (3.7) as a tensor equation. 


—_ (3.12) 


where we have defined the contravariant input vector as 


* 
| 


x(k) x(k—1) --- x(k-A) --- x(k—N)}? (3.13a) 


= xT (3.13b) 


This choice for x has the effect of truncating the series to N+1 terms. The 
symbols H. A) Hy a, --- represent the components of covariant tensors of order 
O. 1, 2, etc.. respectively. 

Examining equation (3.12) we make particular note of the following two 


aspects: 


(1) The dimension of the vector space is related to the memory order of the 
system. The vector x has N+1 components implying that the nonlinear 
system contains no more than N delays. This fact is explicit in the way the 


input vectors have been defined. 


—_— 
lo 


The nonlinearity of the expansion is provided implicitly by the tensor outer 
product operation. It is the outer product of the vector x with itself that 


makes the higher order (2 and up) terms nonlinear. 
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Fig De. : 
igure 3.3: Alternate Second Order Triangular Volterra Kernel Domain 
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These observations provoke speculation about the possibility of an 
alternate formulation where the roles plaved by the dimensionality of the vector 
space and the outer product operation are interchanged. 

2. Alternate Tensor Formulation 
Consider a parametric description of a curve in V®*!, a p+1 dimensional 


vector space, given by; 


x°(k) = [x(k)]" 
x¥(k) = [x(k)]" 
pd eit : (3.14) 


xP(k) = [x(k] 


where the superscript in parenthesis indicates exponentiation. Given any x(k). 
the components x*(k), A = 0,....p define a point in V’*!. As x(k) varies with k, the 
vector x(k) will describe a curve in V®*'. We define the components of another 


vector in V°*! as 

ae | es 0 nee | oe 8 he (3.15) 
Similarly. the N-th such vector can be defined as. 

x’N(k—N) = ix(k-N}O™ (3.16) 


The vectors in equations (3.14). (3.15). and (3.16) will be referred to as 
observation vectors. Although, at this point they only depend on past and 
present system input values. later. in Section D, they will be generalized to 
include past outputs as well. The input and output measurements represent the 
system observations, hence the name. 

Using the theory developed in Chapter 2. concerning multilinear 
functionals, the following mathematical relation 1s proposed as a model of a finite 


order. finite memory nonlinear system: 


§(k) = x°%(k) oO MR-N)H,. (3.17) 


N 
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where H,, ..,, is an (p+l)-order covariant tensor representing a (p+1)-linear 


N 
functional. This covariant tensor performs a role similar to that of a Volterra 
kernel. 

This model is equivalent to the p+1 term Volterra series (equation (3.7) 
or (3.12)), although it does contain many additional terms. It has the advantage 


of notational simplicity. It replaces a p+1 term series with a single term and 


requires the specification of only a single composite kernel, H,....,,. As will be 


shown in Section D, equation (3.17) is considerably more general than equation 
(3.12). It will be shown that Wiener type models can be obtained from equation 
(3.17) by a simple coordinate transformation. Other choices of observation 
vectors yield Autoregressive (AR) or even Autoregressive-Moving Average 
(ARMA) type models. 

In order to illustrate the correspondence of equation (3.17) to the 
standard Volterra type series expansion of equation (3.12) we present a simple 
example. 

3. Example 3.1 

Consider the truncated Volterra series expansion corresponding to 

equation (3.12): 


y(k) = H+ Hx? = Se es BSE (3.18) 


where \1, = 0,1 and A, = 0,1, and 














We may explicitly write out all the terms implied in equation (3.18). This yields 


ik) = 1 
+ Hpx(k) + H,x{(k-— 1) 
+ Hoox(k)x({k}) + Hoix(k)x(k—1) 


+ Hyox(k—1)x(k) + H,yx(k—1)x(k-1) (3.19) 
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Now consider a model of similar order given by equation (3.17). 
A A 
y(k) = Hy yx °(k)x “(k-1) (3.20) 
where Aj, A, = 0,1,2 and 


x°%(k) = [x(k)}° 


x k—1) = [x(k-2)]?” 


The tensor product, x °(k)x"'(k—1), appearing in equation (3.20) results in a 


contravariant tensor of second order. Its elements are 


1 x(k—-1) x()(k—1) 
Ix’(k)x (k-1)} = | x(k)  x(k)x(k-1) — x(k)x(k-1) (2a) 
x(k) x) (k)x(k—1) x!) (k)x()(k—-1) 

We notice that all the elements on or above the southwest- northeast diagonal 
are included in equation (3.19). the terms below this diagonal are not. This new 
form has not discarded any terms present in the latter version and so cannot 
represent a loss of generality. The extra terms that are included in this new form 
do not pose anv significant problem. If the system does not contain terms 
involving these particular elements then the corresponding term in the kernel will 
go to zero during the identification process. Certainly in some classes of systems 


these additional elements will be important. 


C. DISCRETE NONLINEAR SYSTEM IDENTIFICATION 

In Section B a tensor formulation of the discrete nonlinear modelling problem 
was introduced. In this section methods for kernel identification are presented. 
All the methods that are discussed are statistical in nature and utilize a least- 
squares approach of parameter estimation. It will be shown that familiar 
inethods used in linear systems modelling can be extended to handle the 
nonlinear case. In the first section a tensor equivalent of the normal equations. 
which can be solved for the unknown system parameters. is derived. Several 
recursive (in time) solutions to the problem are then presented. Finally. 


simulation results are offered to illustrate the effectiveness of the algorithms. 
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1. Derivation of Normal Equations 
In Section B (equation (3.17)). the following tensor model for nonlinear 


systems was proposed: 


p(k) = x(k) --- x S(kK-N)H,.. (3.22) 


N 


where y(k) represents an estimate of the system output and the x (k-i)’s are 
components of contravariant observation vectors, defined in equations (3.14) 
through (3.16). 

Consider the analysis model of Figure 3.4. This diagram represents a 
conceptualization of the system identification process. The assumption is made 
that the unknown system can be represented by an equation identical to the 


model equation, (3.22), where the system parameters H,,...,, are unknown. The 


parameters of the model are adjusted to best match the actual system 
parameters. A convenient measure of how well the model represents the actual 
system is the mean-square error or MSE. The MSE is a quadratic function of the 
model parameters which implies that there exists a unique minimum. or optimal 
solution, to the problem. Minimizing the MSE yields a linear set of simultaneous 
equations which can be solved for the unknown model parameters. In addition. 
the quadratic nature of the MSE allows steepest descent type. adaptive 
algorithms to be used. | 

The error signal, e(k). defined as the difference between the actual 


nonlinear system output. y(k). and the output of the model. §(k), is given by 
e(k)= y(k)- 9(k} (3.23) 


The mean-square value of this error signal is given by 
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Figure 3.4: Svstem Analysis Model 
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where E{ } indicates expectation. The optimal set of parameters can be found by 
N 


minimizing the MSE, E{e?(k)}. with respect to the parameters H, 4,. The 


gradient of the MSE is formed by differentiating with respect to the unknown 


parameters. 


se{eu} 


OH,... 


0° AN 


= 2 _ x(k) ee x *(k-N)H,, . -eayll—x" °(k) ms <n) (3.25a) 
Setting this equal to zero yields 


ay (Kx"*h).x"(4-N — x"%k)..x S(k-N)x"%(k)...x"8(k-N)H,, | 


=O (Jo) 


Or. 
efelte ae wise). kd) hI, = efvtis"e).x"™k-i} 


[aeeac} 


The expression, (3.25c).,is a tensor equivalent of the Wiener-Hopf or 
normal equations. As will be shown in Section D. these equations can also be 
used to represent the Yule-Walker equations if a different choice is made for the 
observation vectors. x(k-i). The term on the left-hand-side of equation (3.25c) 1s 
a nonlinear extension of the autocorrelation matrix. This contravariant tensor of 
order 2(N+1) includes various higher order correlations as well as the second 
order ones which arise in the linear case. The expectation on the right-hand-side 
represents a cross correlation between the output and various linear and 
nonlinear functions of the input. 


The minimum MSE can be determined by substituting (3.25c) into 


(3.24b). This vields 


ol 


+ Hee Bert) x een stey atv} (3.26a) 


= eb"0} - efi ee Cush, Hye (3.26b) 


Equation (3.25c) represents (p+1)'%*) equations in as many unknowns, 
and so can theoretically be solved for the unknown parameters H,....,,. 
However. in practice the number of computations required to perform the matrix 
inversion becomes unwieldly. An nxn matrix inversion takes on the order of O(n’) 


31N*1)) operations will be required in 


operations |Ref. 27: p. 58], therefore, O(\p+1' 
the solution of (3.25c). In order to make the task manageable. alternate 
algorithms that avoid direct matrix inversion, must be employed to solve the 
normal equations. 
2. The Least Mean Square (LMS} Algorithm 

The Least Mean Square (LMS) algorithm has successfully been emploved 
in the solution of a variety of linear problems ‘Ref. 28]. It is a recursive 
algorithm based on a gradient. steepest descent type of strategy. The inean- 
square error is a hyperparaboloid which is concave upward. Steepest descent 
algorithms strive to descend down this quadratic surface. towards the minimum. 


by making adjustments to the unknown parameters proportional to the gradient. 


This can be expressed matheinatically as 
Hy. | ag (RE) = et) a eee) (3.270) 


where H, ,,(k) is the value of the model paraineters at the k-th time instant 


and pw is a paraineter which controls the convergence of the algorithm. The 


symbol. . (k}. is used to represent the gradient. 


The actual value of the gradient is given. according to equation (3.23a). 


by 


VA, ay(k) = 2E{e(K) x8 x Mk | (ome) 


The LMS algorithm approximates the MSE by the instantaneous value of the 


error squared. The approximate value of the gradient is 


~ 


Ty -++ay(k) = — 2e(k)[x°°(k) «>: x *(k—N)] (3.29) 


where the circumflex indicates an estimated value. Using this approximation in 


equation (3.25) yields 


Hy, ag(K+1) = Hap. a,(k) + 2we(k)x°%(k) «+ x “(k—N) (3.30) 


0 


Equation (3.30) gives a straight forward method of determining the 
system parameters. It involves no matrix inversion nor does it require the 
correlation tensor to be known. These two properties make the calculations 
required at each iteration very simple. so that it is possible to perform them in 
real time. Jt can be shown that the LMS algorithm converges to the optimal 
solution [Ref. 29: p. 578]. 

For the linear case. the algorithm will converge as long as the parameter. 


x, is chosen to satisfy. [Ref. 29: p. 578). 





Or fi =. (3.31) 


CA G4 


where 4,,, 1s the largest eigenvalue of the correlation matrix. Since the 
correlation matrix is positive definite. all its eigenvalues are greater than zero. 
Therefore. the trace of the correlation matrix will always he greater than the 


largest eigenvalue. The following condition will. therefore. cusure stability. 


1 . 
Vg oes 
Tr correlation matrix 


i) 


For the nonlinear case this translates to 


l 


. 7 : aa : ANGEL -\ Aug An, , \ (ean 
Per ee [k) x fe Na “(k)...x (k-N) | 
A y= 0 Ay, = 0 


In practice a value considerably smaller than the maximum permitted by 
equation (3.33) is used to give slower. more accurate convergence. 
3. Recursive Least Squares (RLS) Algorithm 

The recursive least squares. or RLS algorithm, is similar to the well 
known Kalman filter except that it is used to estimate model parameters rather 
than the system state. The tensor form of the normal equations is not suitable 
for RLS implementation. The elements of the correlation tensor must first be 
rearranged in a two dimensional matrix. This can be accomplished by replacing 
the tensor outer products appearing in equation (3.25c) with matrix Kronecker 
products. The Kronecker product of an nx»m matrix, A = ‘aa .> and a pxq 
Matrix. = ib 


yy,/: 1S AN npxmq matrix given by [Ref. 30.31], 


aiBoa.Bi ... anB 
a7Bo ayB -.-- a2,B 

ASB=j . . | (3.34) 
a,iB a,2B ae 


where the symbol * is used to denote the Kronecker product operation. 
In order to rewrite the normal equations. (3.25c). in this matrix form. the 


covariant tensor of system Darametlers, Hy, must be put into a vector form 


An? 
by a lexicographic reordering of the elements. The resulting parameter vector is 


denoted H. The normal equations become 


ef xt (oo et = Ni sc Re ear ee ice yn 


chin KU ee Scie xy (Sie 


If an assumption of crgodieity is made then the statistical averages in (3.35) can 


‘ 
< 


be replaced with time averages. 


1 tro ce a ae 
Xiki kok) | eel {meth 
7 au (kK) (k] lim = dy )X(h} 


jim ay 
=e 











where X(k) = x(k) <-°-° © x(k N). For computational purposes equation (3.36) 


is approximnated as 
K-1 - K-1 
dy X(k)X7(k) = SY) y(k)X(k) (3.37) 
k=0 k=0 


As a further notational convenience, the following two definitions are made, 


X*(0) 
X*(1) 
——| (3138) 
*(K) 
and, 
y (0) 
Bil.) 
y-| | (3.39) 
y(K) 


Using these definitions the normal equations can be written in the 


compact form 
XJX,H = XY (3.40) 


This last from of the norinal] equations is precisely the same as the one used by 
Goodwin and Payne 'Ref. 32: p. 176) for the derivation of the RLS algorithm. 
The derivation involves the use of the matrix inversion lemma /Ref. 33: p. 
247| which replaces a matrix inversion at cach iterative step by a simple scalar 
division. It is this simplification which yields the efficiency of the RLS algorithm. 


The RLS equations are |Ref. 32: pp. 176-177). 


Hy, ie Ore (y Uh) XI (K - 1) Hy) (3.41a) 


Pe XK - 1) 


a (aed XK 1) 


(3.416) 


Py, = (1-0, eee (3.41c) 


=X) xe (3.414) 


Equation (3.4la) deserves some comment. The term in the parentheses 
is norma!ly known as the innovation. It is a scalar which represents the new 
information gained in the latest measurement. The term X1(K+1)Hy is an 
estimate of the current output given the new input measurement and the old 
estimate of the system parameters. The innovation is thus an error signal. The 
term in front of the bracket, Qx,,, is a gain vector. It gives an indication of how 
much faith is being put in the new information. If the gain is small. then the 
new estimate will be essentially the old estimate. Conversely, if the gain is large 
then the new estimate will depend to a great degree on the new information. 
When the algorithm has converged. the gain will be close to zero. If the system 
parameters change for any reason. the algorithm will not detect the change since 
it is ignoring the new-information. In order to circumvent this problem a 
weighting can be applied to the data. so that new data is artificially emphasized. 
Exponential and rectangular windows have been successfully emnploved for this 
purpose. The time constant used in the case of the exponential window is often 
called a forgetting factor since it ensures that the algorithm's memory is finite. 
Use of windows will not be discussed further in this dissertation. The interested 
reader should consult reference [Ref. 32: pp. 179-185]. 

4. Simulation Results 

In order to investigate the validity of the theoretical results. several 
colulputer simulations were performed. Two examples are presented. representing 
two different classes of systems. Many other systems were also tested. however. 
the results presented are typical of those obtained from all the simulations. The 
FORTRAN programs written allowed nonlinearities of up to fourth order. but 
were limited to systems lnvolving only zero or one delay. 

The first example was chosen to correspond exactly to the model 


equation (3.17). The system was excited by white, zero mean. uniform noise. 
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The problem was to identify the system parameters from only input and output 


measurements. The unknown covariant H,,,, tensor was chosen to be 


2 —4 03 —.7 0.0 

Poe TOOws oA 9 0.0 
[Fa.a,| =o =993 .7. 010 (3.42) 

43 .81 -.05 .4 0.0 

Oo 007 016 0.0 0.0 
Therefore, the output equation consists of 16 nonzero terms, ranging from 
2x (k)x)(k-1) up to .4x')(k)x®(k-1). This model will be called System I in the 

sequel. 

The other type of nonlinearity tested was one that was known to require 
an infinite number of terms. The particular example chosen for this was the unit 
step function which simulates a saturation type of nonlinearity. It is a convenient 


choice since an analytical solution can be calculated. 


The unit step function has different H,,  ,, tensors depending on the 


order of the model chosen. This is a result of the chosen coordinates not being 
orthogonal (this will be clarified in Section D.) The second. third and fourth 


order {nonlinearity order) models are given bv 


aes: 
|Ha,| = p 7 J 





ls Pavezs 0.0 - 1.09875] (3.45) 








| i oe 35 
[i | | 2 9 | 


= E 1.40625 0.0 —1.09375 0.0] (3.43c) 
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Notice that the unit step function has no memory. The expansion contains only 
odd powers of x{k) and a constant term since the unit step function can be 
written as a constant plus an odd function (the signum function.) We will refer 
to these three models collectively as System II in the sequel. 

The direct solution of the normal equations was tested on these two 
models as was the LMS algorithm. To date the RLS algorithm has not been 
verified. 

a. Direct Solution of Normal Equations 

To verify the direct (matrix inversion) method of solution of the 
normal equations, a FORTRAN program was written which estimated the 
correlation tensors, and performed the required matrix inversion. The program 
Was written so as to allow the number of points used to approximate the 
correlations to be varied. The maximum power of x(k) was also made to be 
adjustable so that the effect of over or under modelling could be studied. The 
final adjustable parameter was the magnitude of the uniform, zero mean. white 
noise that was used to excite the system. Adjusting this last parameter affects the 
range over which the resultant model is valid. In an actual application. 
something about the range of expected system inputs would have to be known in 
order to select a good value for this parameter. 

The results for System I are given in Table 3.1 for several 
combinations of the three variable parameters. The results are remarkably good 
even with as few as 30 input samples. Since no measurement noise was 
introduced this is perhaps not surprising. 

Overmodelling did not present any problems. The additional terms 
were identified by the algorithm to be essentially zero. This is evident in Table 
3.lc. Undermodelling did introduce some inaccuracy. The coefficients identified 
are significantly different from the ones in the unknown system. However. the 
coefficients identified should represent the best second order approximation to the 
system. which will not in general be the same as the second order coefficients 


contained in the third order model. 





a. 30 data points were used 
The maximum power of x(k) was 3 
The noise was uniform on (-1,1) 


200 
500 
Ayal = | o10 
430 


b. 500 data points were used 
The maximum power of x(k) was 3 
The noise was uniform on (-10,10) 


19993 —.39952 


} 49990 
Aaa, ~ | 010006 
43000 


c. 500 data points were used 
The maximum power of x(k) was 4 
The noise was uniform on (-1,1) 


. 20000 — .40000 
.90000 .39000 
IHW, = .010000 1.30000 
.43000 .81000 


—.48140E-06 —.43483E—06 


d. 500 data points were used 
The maximum power of x(k) was 2 
The noise was uniform on (-1,1) 


—.400  .030 
.350 ~=.110 
1.300 —.330 
.810 —.050 


20616 ~ .80030 
‘H,,:= | .75931 1.5108 
~.033444 1.6619 


Omit 


.343539E—05 


700 
900 
700 
400 


.030016 —.70001 


.380014 .11001 .90000 
1.3000  —.383000 .70000 
81000  —.050000  .40000 
.O3000 — .70000 
.11000 .90000 

— .383000 . T0000 

— .050000 .40000 


.52571E—06 


~ 26356E-06 
~ 10016E—06 
33811E—05 
— .83230E—07 
~.43294E-05 


TABLE 3.1: System I Results Using Full Matrix Inversion. 
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The results for System IJ are shown in Table 3.2. Significantly more 
points were required to accurately identify System IJ than were needed for 
System I. The unit step function cannot be represented exactly by a finite 
number of polynomials so it is not surprising that the solution is not precise. 
There is another factor that contributes to the poor accuracy. The choice of 
functions used as coordinates, equation (3.14), leads to ill conditioned equations 
(see Strang [Ref. 34: p. 135].) This problem will be corrected in Section D. 

b. Simulation Using LMS Algorithm 

To verify the LMS algorithm a FORTRAN program was written, 
which used equation (3.30), to adaptively identify the covariant tensor Hae. 
The program allowed the convergence parameter. », and the magnitude of the 
uniform. white noise to be varied. The results of the simulations are presented in 
Table 3.3 and 3.4 for System I and II respectively. Equation (3.33) was used to 
bound the convergence parameter. ». The input excitation noise was chosen to 
be uniform on (-1.1) resulting in a bound for the convergence parameter of 
On 0 ooo. 

In general convergence was slow. The linear and "close to linear" 
terms were identified most rapidly. The highly nonlinear terms. involving high 
powers of x(k) or x(k-1) and their cross terms. were last to be identified. Their 
accuracy never reached that of the lower order terms. The algorithm was very 
sensitive to the setting of the convergence parameter. ». A value smaller than 


the bound predicted above was used to achieve satisfactory perforinauce. 


D. GENERALIZED COCR DIN 3A be = Ub 

In Section B. a coordinate system was introduced that was closely related to 
the Volterra series (equation (3.14).) This system was subsequently used in the 
remainder of Section B and in Section C. There is really litth: motivation for 
choosing this particular set of coordinates. In fact there are very compelling 
reasons to search for other sets of coordinates. The set (3.14) can lead to poorly 
conditioned sets of equations (see Strang Ref. 34: p. 135]). a fact that was 


mentioned in the last section. In higher order systems this can become a serious 


6O 


a) 


. 1000 data points were used. 


The maximum power of x(k) was 3. 
The input noise was uniform on (-1,1) 


49278 1.4301 ‘0005139 —1.be39 
WI ‘ —.037083 —-.10437 .067337 .18706 
eal (029155 —-092552 = 7085273 .10310 
084649 .080036 —.13334 —.15251 


. 15000 data points were used. 


The maximum power of x(k) was 3. 
The input noise was uniform on (-1,1). 


50380 1.4131 —.0093646 —1.1029 

~ 0064558 —.024451 .030395  .026762 

Haas) = | — 0016623 —.021414 .0044444 030779 
00015024 .024709 —.029871 —.018542 


. 500 data points were used. 


The maximum power of x(k) was 2. 
The input noise was uniform on (-1,1). 
47416 Tti4G 00062628 
H,,. = | 0049366 —.0018654 = 0025596 
auz02 199) — . 10098 07699] 


. 5000 data points were used. 


The maximum power of x(k) was 2. 
The input noise was uniform on (-1.1). 


= 


A9480 i oaae G2 (lac 
ya = | 017591 - 001427 O24aTi4 
@5590 0019598 (eae 


TABLE 3.2: system I] Results Using Full Matrix Inversion. 
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The maximum power of x(k) was 3. 
The noise was uniform on (-1.1). 
The convergence factor, 7. was chosen to be .2. 


a. After 100 Iterations 


21191 —.54901 .23961 —.58084 
49735 .55747 -—.22910 .62263 
Hag] = | 11778 1.3872 —.42404 85335 
43573 41638 —.29008 .45359 


b. After 300 Iterations 


19133 —.45457 —.089402 —.69479 
52472 .65421 .055704 = .71551 
032839 1.24838 —.34090 .75055 
44209 .55546 —.086633 .53219 


i =a 
‘Hy a,] = 


c. After 500 Iterations 


.18995 —.35146 088102 —.68228 
01438 .55292 047074 66728 
0283802 1.2505 —- .28989 


Aya] = 799 
MoO0T mE OO2co OS2106 50355 


d. After 1000 Iterations 


20000 39401 0289125 - 100.) 
.00360 =. 46156 1090] 14664 
017045 1.2743 —.32721) 75052 
15879 64126 = (O10G10 {Gia 


Le), 


e. After 1800 Iterations 
19456 -—.38962 .030956 7 
48923 eA ei iy 4 HUG 
AgAy O07 90a2 122786 — 939206 i 


foc 67660 Oaoae8t 61097 


TABIBS.3: System 1 Results sine DNs Vicor 
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The Maximum power of x(k) was 3. 
The input noise was uniform on (-1,1). 
The convergence factor, », was chosen to be .15. 


a. After 100 Iterations 


27229 64437  —.029974 —.30266 
— 021558 —.13134  —.16022 —.10782 
(Hy... = | 962328 049563 —.069966  —.17809 


.26624 —.00231392 —.093028 —.048308 


b. After 300 Iterations 


49886 .96021 —.11685 —.88415 

~.010596 .10195 —.060408 —.12604 

Hag = 1 036964 15373 -.17940 —.32627 
14847 059296 —.049165 —.10986 


c. After 500 Iterations 


54880 1.3925  .010187 —.88497 
__ [-.091216 .14549 10978 ~ .0080069 
Anos = | 981292 .23020 —.14849  — 3387] 


0091398 .021484 -.012741  -.071461 


d. After 1000 Iterations 


mo2o0 I.o2208 9=.034331 - 1.0004 
7 ; .13698 MGUZ70 leods O36125 
eA 0AY 16666 .16941 —-.096854 —.35439 
068181 018017 084221 (0062654 


e. After 1700 Iterations 


“eos? 1.3719 —.058312 ae95098 
PoosteO25o2s O12228) 027707 

ome | 096072 .20474 —.0R4G94 —.2153 
Mesos 047110 =.093643 025036 


H, 


TABLE 3.4: System II Results Using LMS Algorithm 


problem. It is shown in this section that proper choices of coordinates lead to 
diagonal correlation matrices so that no costly matrix inversions are required to 
solve the normal equations. This makes the identification process almost trivial. 
It was stated in Section B.3 that the input observation vectors could be 
considered to describe a curve in a hypothetical (p+1)-dimensional space. The 
system output ‘is estimated based on this curve (equation (3.17).) This chapter 
generalizes this idea to include cases where we have two curves, one depending on 
present and past inputs, the other depending on past outputs. Finally, the 
chapter concludes with several non-trivial numerical simulations. 
1. Choices of Coordinate Systems 

A point in (p+l1)-dimensional space is defined by the variables 
x°, x? ..., x®. In Section B.3 (equation (3.14)) these variables were chosen to be 
parametric functions of the variable x(k), the system input. The particular 
choice presented there was chosen to correspond to the Volterra series. It is 
desirable to pick a system of coordinates that ensure a diagonal correlation 
tensor. as this allows solution of the normal equations without matrix inversion. 


A tensor T is diagonal if its components obey the rule 


pio aw 2 2 2 (3.44) 


otherwise 
0 


> {Ores A «eee ee abe Ae =X. 
Ky 0 ea 1 Nel chia N 


Note that only tensors of an even order (possessing an even number of indices) 
can possibly be diagonal. In general components satisfying the upper condition of 
equation (3.44) are called diagonal elements, or diagonal components. 
Components that are not diagonal are called off diagonal. 

Two conditions are required in order to ensure a diagonal correlation 
tensor. The first is a result of the following theorern. 

a. Theorem 3.1 

If a set of functions { f(x). f(x). .... fy(x)} are defined with the property 


that 
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Ky for A = 
ne x}w(x}d <= {% ee a ‘ (3.45) 


where w(x) is a positive weighting function, then the set of random variables. Z*. 


defined as 


Z° = fo(X) 
Z' = f,(X) 
(3.46) 
Z? = £,(X) 
where X is a random variable with probability distribution 
px(x) = Cw(x) (3.47) 


are uncorrelated. In equation (3.47) the constant. C. and weighting function. 
w(x), must be chosen so that px(x) satisfies the definition of a probability density 
function. 


b. Proof 


e(va'| = e{0(%) (3.48a) 


ox 


= f S096, Ce)px(x)ds (S.48b) 
Joona x)Cw(x)dx (3.48c) 
- C f f,(x)f,(x}w(x}da (3 48d) 


Substituting equation (3.45). yields 


ie fora 
ae ie fOr 
cfr Q ford 


Choosing a set of functions in accordance with equation (3.45) and 


BP (onto) 
pe 


td 


ensuring that f)(x) is a constant function (ie: is a constant) is the first condition 
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that must be satisfied in order to obtain a diagonal! correlation tensor. The 
system must be excited with samples drawn from a random process with 
probability distribution given by (3.47). In addition. if the random process is 
chosen to be zero mean and strictly white (ie, independent which implies 
uncorrelated), then the correlation tensor will be diagonal. This last condition is 


equivalent to requiring that 


E4x(m)x(n)? = 6(m—n) (3.50a) 

p(x(n),x(m)) = p(x(n))p(x(m))} (3.50b) 
and 

Eyx(m); 0 (3.50c) 


where x(m) and x(n) are two input samples taken at times m and n respectively. 
It is a straight forward matter to show that if these conditions are met. the 
correlation tensor will be diagonal. 

The conditions presented above imply that different sets of 
coordinate functions should be used depending on the probability distribution of 
the noise used to excite the system. Different noise distributions have the effect 
of weighting the error differently. Consider that a Gaussian noise will contain 
samples of all amplitudes while uniform noise is bounded. If. for example. the 
system contains a saturation type of nonlinearity. the uniform: noise may not 
detect its presence if its maximum amplitude as not sufficiently large. 
Theoretically. Gaussian noise contains samples of all amnplitudes and will excite 
all modes. On the other hand it may he known that tlic system: imput never 
exceeds a certain maximum value aud so « bounded input will be suitable. 

If two models of a svstem are constructed. in two different coordinate 
syetetus but using the same input noise (only one set can possibly lead to a 
diagonal correlation tensor) then the two solutions will be equivalent. One 


solution can be transformed into the other by performing a change of coordinates 
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in accordance with equation (2.22). Therefore. it always makes sense to identify 
the system using the coordinate systern which leads to a diagonal tensor. Other 
representations can then be calculated as desired. Solutions obtained by using 
different input noise density functions are not equivalent even if the coordinate 
systems used were the same. The effect of the different distributions is to weight 
the errors differently. A uniform input noise will weight the errors equally, while 
a Gaussian noise will emphasize the importance of errors made for small inputs. 
In general, transformations between solutions obtained using different excitation 
noise distributions cannot be found. Choosing an appropriate distribution 
requires knowledge of the expected system input signals. 

The Hermite polynomials lead to a diagonal correlation tensor if the 
input is white. zero mean, Gaussian noise. Similarly, Legendre polynomials 
should be used in the case of uniform noise. It is convenient to normalize the 
coordinate functions so that the diagonal components of the correlation tensor are 
all ones. To identify the system parameters. only the cross-correlations of the 
right-hand-side need to be calculated. This fact was first discovered by Lee and 
Schetzen [Ref. 10}. 

2. Recursive Models 

Recursive models have been used very successfully for modelling lnear 
systems [Ref. 35]. Among their advantages is infinite memory. and the ability to 
model a system without knowledge of its input. The latter property allows these 
models to be em:ploved in such areas as speech processing where input signals are 
difficult or impossible to measure. The assumption is made that the input 1s 
white noise. 

Recursive discrete-time nonlinear expansions have also been proposed 
(see for example /Ref. 14.15.16]). Recursive models posses infinite memory. and so 
mlay require fewer terms to accurately represent long memory systels. However. 
nonlinear recursive models also posses infinite nonlinearity. To understand why 


this is true. consider the following example. 
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a. Example 3.2 


Consider the following recursive, discrete. nonlinear system 
y(k) = ay?(k—1) + u(k) (gee 


where u(k) = bé(k) (where 6(k) is a unit sample) and where a and b are arbitrary 


constants. 
The output (y(k)) of the system for k = 1, ..., K is 
AD (3.52a) 
y(1) = ab? (3.3215) 
y(2) = a(ab*)’ = a®b* (3.52c)} 
CS) = cP (3.52d) 


It is clear that the nonlinearity of the system increases with time. Unlike a linear 
system. stability in a recursive nonlinear system is determined not only by the 
svstem parameter. a. but also by the input function. It is also difficult. in 
general, to predict for what classes of input a particular system will be stable. 
By analogy to the linear case we will refer to these types of models as Auto- 
Regressive or AR. 

It is possible to also expand the output as a combination of both past 
and present inputs and past outputs. This type of model was proposed by Parker 
and Perry ‘Ref. 13]. We will refer to this type of model as an Auto-Regressive. 
Moving-Average or ARMA model. It will have the same stability problems as 
does the AR model. 

Equation (3.17) can be used to model an AR nonlinear systems if the 
proper choice is made for the observation vectors. Using the same coordinate 
functions as given in equation (3.14) but using y(k - i), i = 1,....N, as an input 


parameter, yields appropriate observation vectors. That is. 
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y (ki) = 'y(k- yp (3.53) 


defines the components of the i-th observation vector. The model equation. 


equivalent to equation (3.17), becomes 
d d - 
mie = y (k~1)...y “(k-N)H, ...,, (3.54) 


where H,,...,, is again a (p+1)-th order covariant tensor containing the system 


parameters. This model is an extension of the familiar linear autoregressive, or 
AR model. 
The normal equations derived for the nonrecursive case can be used 


to solve for the H,,. ,, tensor in this case as well. The identification process is 


described in Figure 3.5. The assumption is made that the system is recursive and 


driven by a white noise, u(k). Its output can then be described by 
rx o = 
y(k) = y(k-1)..y “(k-N)Hy,..- a, + u(k) (3.55) 


This output signal is delayed and fed into the analysis model. The error signal 


e(k) is given by 


e(k) = y(k) — §(k) (3.56a) 


y k-tj..y “(k-NJH, ©, ~ u(k) — y (k- 1)... “k—-N)H, 


AN 


(3.56b) 


When the model parameters exactly equal the actual system 
parameters the error signal will equal the input white noise. For this reason the 
analysis model is often called a whitening or bleaching filter. The analysis 
model is nonrecursive. It uses past values of the system output to make a 
prediction of the present svstem output. The normal equations (3.25¢} apply to 
this situation with the observation vectors. x(k-i). replaced by the vectors defined 
in equation (3.53). The normal equations for the recursive nonlinear model can 


be written as 
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Figure 3.5: Recursive Model Identification 
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phy th) 9 UN9 TRI) "ANI, dy 


= efyhiy"k=1)-9"%i-] (3°57) 


These equations are a tensor equivalent of the Yule-Walker equations. The 


corresponding mean-square error 1s 


e(e00} = eb") = Pp (ks Le vanaf, . (3.58) 


In this case it is not obvious that one set of coordinates will yield 
better results than another. The correlation tensor appearing on the left-hand- 
side of equation (3.57) cannot be guaranteed to be diagonal since the probability 
distribution of y(k) cannot be controlled. Techniques must be devised which 
choose the coordinate system "on line" as the statistics of y(k) are determined. 
For example. in the linear lattice the coordinate system is chosen by 
orthogonalizing the sequence y(k) using a Gram-Schmidt procedure. In this way 
the coordinate system is not determined until run time. Extension of these ideas 
is left until Chapter 4. 

The tensor model presented is also suitable to represent a nonlinear 
ARMA model. This type of model takes into account all available information 
(input and output) and so should be more accurate. It may also lead. in some 
cases. to a lower order solution than either the AR or MA model. An ARMA 


tensor model can be written as 
A A a 2 
i= x “(k)..x “(k-M)y (k-1)...y “(K-N)Hy,. aye, ny (3.59) 


Relations analogous to (3.57) and (3.58) for the normal equations 
and the minimum mean-square error for the ARMA model can be obtained in a 
straight forward manner. 
3. Simulation Results 
The concepts developed in Section D of this chapter were verified using 


FORTRAN simulations. The coordinate functions were chosen to be the 
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normalized Legendre polynomials because of the ease in generating uniform noise 
on a digital computer. By assuming that the correlation matrix was diagonal. 
the solution to the normal equations was determined without matrix inversion. 
This solution was verified by generating the correlation on the left-hand-side of 
the normal equations and performing the required matrix inversion. The LMS 
adaptive algorithm was also tested using the Legendre polynomials. 


The Legendre polynomials used are given by 


Aes eile (3.60a) 
f(x) = Vax (3.60b) 
ba) = ~ (| (3.60c) 


Hea) / == (x(3) =x) (3.604) 


These functions are normalized so that the correlation tensor will have ones along 
the diagonal when excited with zero mean white noise which is uniform on (-1.1). 

The simulation models used were similar to those used in Section C. The 
coefficients for System I were identical to those given in (3.42), although thi- 
time they are coefficients of Legendre polynomials so they do not represent the 
same system. Svstem I] was again a unit step function which has a tensor 


representation, in terms of the Legendre polynomials given by |Ref. 26: p. 1526}, 


H, age yee (3.61a) 
2 4 80 


los GeesOls ny 0.165359] (S.61b) 


The results of the simulations for System I] are presented in Tables 3.5 
and 3.6. Results for System I] are given in Tables 3.7 and 3.8. In all cases it is 


obvious that the simulation results are very good. 


15000 sample data points were used. 
The maximum power of x(k) was 3. 


a. Solution Without Using Matrix Inversion 


19158 —.38475 .0041187 —.67056 


49067 .38457 .089293 .93948 
Easels =9007253 480995 —.35452 .71923 


.41049 85202 —.074247 .43011 


b. Solution Using Matrix Inversion 


.20000 —.40000 .030000 —.70000 
.90000 = .35000 .11000 .90000 
01000 1.3000 —.33000 .70000 
4800 .81000 —.050000 .40000 


ie) = 


TABLE 3.5: System I Model Parameters Obtained 
a. Using No Matrix Inversion, and 
b. Using Matrix Inversion. 
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The maximum power of x(k) was 3. 
The convergence factor was chosen to be .01. 


a. After 200 Iterations 


14678 —.45852 —.036398 
42720 .25876 .018320 
(Aaa, ~ 1 — 077103 1.1898 —.45388 
34264 .68052 —.16955 
b. After 500 Iterations 
20013 —.40047 .030509 
| | 50038  .35077 .11136 
Aaa, = | 0099732 1.2996 —.32958 
43008  .81014 —.050063 
c. After 1000 Iterations 
20000 -.40000 .030007 
| 50001 .35000 .11001 
ee, ie 010008 2000 32999 
43001 .81000 —-.049986 


bP FA TT 
.83079 
.61647 
.30546 


—.70066 
.90025 
.70026 
99875 


— 0000 
veel, 
[oes 
.40000 


TABLE 3.6: System I Model Parameters Using LMS Algorithma. 
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a. 500 data points were used. 
The maximum power of x(k) was 3. 


No Matrix Inversion Solution 


.52104 .43930 0090070 —.17515 

—.011098 —.0046220 .014464 .020122 

Hag, ~ | .0011160 —.010606 —.010466 .021310 
.040052 —.0071009 .058263 .052566 


Using Matrix Inversion 


00157 .43054 = JIiso, —.16530 
| — .0084508 .0032387 —- .Q000970 —.0006363 
Aaa z= .0017041 —.0009167 .0062532 —.0059126 


— .0066652 —.0024676 .012192 .0014058 


b. 15000 data point were used. 
The maximum power of x(k) was 3. 


No Matrix Inversion Solution 


.50050 43363 —.0043279 —.16772 
; - 0015150 .0003272 -—.0012861 —.000373 
Fy ay) ~ | —.0071103 ~—.0004279 .012933 014358 


nO022021 008209 —.U021332 — 70060074 


Using Matrix Inversion 


SOORYG 43306 — QOO1744 WiGOG 22 
018212 —.0005741 - 0014410 — 0000789 
Hie A036401 .0008143  .0043811 — .0008015 
QO0038990 .0001746 - .00032974 .O0O0EC4I67 


1 3b Tae 
a. Using No Matrix Inversion. and 
b. Using Matrix Inversion. 
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svstem I] Model Parameters Obtained 


The maximum power of x(k) was 3. 
The convergence factor, », was chosen to be .005. 


a. After 300 Iterations 


AG464. «41356 © .0066999_ —.16385 

011790  .012146 .0067764 —.0093346 
[Hawa] = | — 010043 .014180 020212 —.018152 

0014149 —.0065302 .0093879 —.0047952 


b. After 500 Iterations 


50556 43020 —.00098824 —.16796 
a ~.015344 —.010684  .020836 .0098449 
Aya = | 9082813 —.0052115 .0090727. —.010251 


012877 = —.0094530 —.0046525 .0072316 


c. After 1000 Iterations 


49868  .44629 0053940 —.17370 
| 020577 .0091397 —.010215  .0065185 
Haas, = | 019857 020072 ~.0064996 —.010541 


—-0020197 .0031219 — (0064275 3 0027119 


c. After 1800 Iterations 


.00115 AQT27 — 00599 File 16505 
000129 —.006G172 — 012712. J0lo2is 
ie. O065666 .0084999 —.0068403 - 010135 
00735684  .0058007 ell 5s 118235 


TABLE 3.8: System II Model Parameters Using LM{IS Algorithm 
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Py eee veOmemeA LT TICE FILTER STRUCTURES 


This chapter reviews lattice filter theory. These filters were first proposed in 
connection with the linear prediction of speech by Itakura and Saito [Ref. 36] in 
1971. They developed a new approach based on a partial correlation (PARCOR) 
coefficient. Since that time the filter structure they proposed has come to be 
known as a Lattice (or sometimes also called Ladder) filter. The properties of the 
filter have been exhaustively studied by many researchers [Ref. 37]. Lattice 
filters have been successfully applied to problems in various disciplines. 

The great interest in the lattice approach stems from it’s property of 
orthogonality. This property allows the filter to be updated in order, without 
recalculation of all the previous, lower order, filter coefficients. Orthogonality also 
leads to a nice physical structure, a cascade of first order sections, and so is 
appropriate for efficient hardware or software implementations. Finally. it has 
been shown that the lattice owes its robust numerical henecelent: to this property 
of orthogonality [Ref. 38: pp. 128-136]. 

This chapter begins with a derivation of the one-dimensional (1-D) lattice. 
The approach en in the derivation is somewhat novel in that it begins by 
expressing the linear prediction in terms of an uncorrelated error sequence. This 
is regarded as a change of coordinate systems and used to develop the Levinson 
algorithm and the lattice filter. In the following section generalized forms of the 
Levinson algorithm and lattice are derived. It is shown that these also 
correspond to coordinate transformation. Next. the Schur algorithm. which is a 
method for generating the required filter coefficients (Lattice parameters) 
directly. given only a knowledge of the correlation matrix. is reviewed. In the 
following chapter the generalized lattice filter is used to develop new. 
multidimensional (specifically 2-D) lattice filters. In Chapter 6. a new nonlinear 


lattice filter, based upon this generalized lattice formulation. is presented. 
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A. 1-D LINEAR AUTOREGRESSIVE LATTICE FILTER 
We define the N-th order forward error sequence of an autoregressive model 


as 
e®(k) = y(k) — by hy y(k — A) (4.1) 


This can equivalently be written in tensor notation as, 


eV (ke) eahyaa (4.2) 
where 

ly? |? = (yk (kj yt kek) (4.3) 
and 

lhy))} = [1 —h,* —h’ --- —hy’ 0...0) (4.4a) 
where for convenience. we have made all vectors of length K+1 (ie; A = 0. ..., K), 


where K is some arbitrary maximum length. Note that, 
ho. =f : (4.4b) 


The y* can be considered to comprise a coordinate system in an K+] 
dimensional space. The vector y*, represents a single realization from a random 
process. As mentioned in the Chapter I], there will in general be many such 
vectors corresponding to all the possible realizations of the random process. 

Because the error. e*(k). is a scalar (invariant) it must remain unaltered 


regardless of the choice of coordinate system. We may. therefore. write 


e* (eee aces (4.5) 


where the y* are defined as 
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y(k) 
(k—1) 
y(k—2) — hgy(k-1) 


y(k-3) my hyy(k—2)hgy(k—1) (4.6) 


SSIS) hy-zy(k—-K+1) =) hg 'y(k—-1) 
and 
[Ky] = [1 —K, —K, aed —Kyn 0...0] (4.7) 


where the h’s are chosen so that the components, y* . are uncorrelated. The 
stochastic form of the Gram-Schmidt procedure can be used. A discussion of this 


method can be found in [Ref. 39: pp. 382-383]. By uncorrelated we mean 


SY fora = py 
dea = 4.8 
ef ; , otherwise ( 


The reader familiar with lattice structures will recognize the components of y* as 
the backward prediction errors. That is they are the errors in predicting y(k-N) 
from the next N-1 values: v(k-N+1), ..., y(k-1). 

It is a straight forward matter to solve the prediction problem given by 
equation (4.5) (ie; solve for the K’s) because of the uncorrelatedness of the chosen 
coordinate system. Using an approach similar to that presented in Chapter III. 
Section C. a set of normal equations can be formulated for this problem. In this 
case. however. the correlation matrix is diagonal so that there is no inversion 
necessary to obtain the solution. Minimizing the mean square value of the error. 


e*(k). given by (4.5). with respect to the K}\’s. vields 
G--: O (4.9) 


Having obtained a solution to the orthogonal problem. we wonder if it cannot 


be employed to advantage to simplify the calculations required to solve the 
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original autoregressive problem. From equation (2.13). we know that the 


relationship between the old and new components can be expressed as 











ov 
ee = a y* (4.10) 
where 
1 0 0 0 0 
0 I 0 0 0 
Oho 1 0 O 
Ceehe  =he 0 0 
dy? 
SG Prey ; : ; (4.11) 
ee he oe 2h, ih i 
Ooh — he =i 
Now. since 
e*(k) = y*hs = y* KS (4.12a) 
, By" > N 
are i (4.12b) 
then 
ro Soe ae (4 13) 
N= A an -LwvJ 


Equation (4.13) gives the relationship between the autoregressive model 
parameters and the K;\’s. This result wil] be used in the proof of Levinson 's 
algorithm. 

1. Levinson-Durbin Algorithm 

In 1947 in his classic paper [Ref. 40) Levinson developed a inethod for 
recursively solving the normal equations. Beginning with @ zero order solution 
successively higher order solutions are calculated. This algorithm can be used to 
exploit the Toephtz nature of correlation matrices of stationary randoin processes 
in order to reduce the required number of computations. The algorithm as 


presented in this work is actually a simplified version of Levinson's original from. 
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It assumes that the equations being solved are the Yule-Walker equations (|Ref. 
41], see also Chapter 3. Section D) so that the right-hand-side of the equations 
contains only terms that will be present in the correlation matrix (left-hand-side) 
of the next higher order model. This simplification is due to Durbin [Ref. 42]. 
We will refer to this algorithm often simply as the Levinson algorithm, however, 
Durbin’s contribution is acknowledged. 
a. Theorem 4.1 Levinson’s Algorithm (Durbin’s Form) 
The autoregressive model parameters may be calculated recursively 


from the relation 


io) = (hf| — Ky.,Sthy" (4.14a) 
where 

fee oh? bh --- -ho, 10--- oO (4.14b) 
and 

meee 0 -h, —-h, --- —hy, 1 0 --- oO (41Ac) 


The operator S has the effect of shifting the components. h; one position to the 


right. Note that 
hy = 1 
For stationary processes the following simplification applies 
h; =h<,., fordA=0.....N eleld) 


b. Proof (by Induction) 


Using equation (4.13) for the first order model we have 


hy Pewee sa) (4.15a} 
a SG 6) Bea) (4.15b) 
10 -:: 0 K,010--- 0 (4.15¢) 
Sie aK, Silay: (4.15d) 
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so that equation (4.14a) holds for 


N-th order model 
Ih) = {hX— "| — KySihy 
From equation (4.13) we have 


yas 


N_ iN 7 De 
nN = KNEE AA’ = 0, K 
Therefore, 
1 
meee Khe + Khe eee ene 
Gee Koh; 2 Ki hy + oe, 
hy = | ~Ky.1+ Kyhy—2 
—Ky 
0 
0 
Now, for the (N+1)-st order model 


] 
4 r 7 me Lae 
. Is, ny Kens Pica Kho ae Ky,hg 
r a sy re. ; Yl a 3 r a N 
ar kK, eke K3h, Tate K,hy “alectvere aati Ky,h, 


\ 4 
eae 


Nee 


We assume the recursion holds for the 


(4.16) 


(4.17) 


1 


2 


(4.18) 


(4.19a) 
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7 


a 0 
1 : ee 
IK ris Koho = Ksh¢ 9 ST Kana. 4 
i j { N- ~h 
—K, a Kenn is K,h; F665 5h Kyh, ! 1 
7 K os Kn ne (4.19b) 
Sis 
0 1 
0 
0 
0 
\ 
= ‘ht — KyerSlby (4.19c) 


And so the desired result has been confirmed. 

The proof presented does not require that the process be stationary 
and so the condition of (4.14c) is not necessary. One is then faced with the 
problem of how to determine the h}\’s. This question is answered in the next 
section where a generalized form of the algorithm is derived. We note that the 


condition of (4.14c) implies the following 
Me Ke = he | (4.20) 


or that the second column of the matrix given in equation (4.11) contains all the 


lattice coefficients. 
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The following significant observation about equation (4.14a) is made. 
The transpose of the term within the second bracket on the right hand side 
corresponds exactly with the rows of the coordinate transformation matrix of 
equation (4.11). The Levinson algorithm, therefore, represents a recursive 


method for finding an orthogonalizing coordinate transformation. 


2. The 1-D Lattice Structure 

It has been shown, in the previous section, that the autoregressive model 
parameters can be calculated in a recursive manner using the Levinson algorithm. 
It was also shown that a model could be built in an orthogonal coordinate system 
where the model parameters were given by the K,-’s. The question arises, can a 
filter structure be devised which represents y(k) in the orthogonal coordinate 
system? The answer is affirmative’ It will be shown in this section that the 
Lattice filter is the required structure for the stationary case. A more general 
solution is presented in Section B of this chapter. 

The desired result is obtained in a straight forward manner by 
multiplying both sides of equation (4.14a) by y*. This yields 
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eN**(k) = e*(k) — Kyi rX(k-N-1) (4.21) 


where the quantity r‘(k—N-1) = y*? evaluated at A’ = N+1. As was mentioned 
earlier in this chapter. the quantity r‘(k-N-1). is generally known as the 
backward prediction error since it corresponds to the prediction of the point y(k- 
N-1) from the N future points y({k-1). ....v(k-N). 


Assuming stationarity. we can use (4.14c) and (4.14a) to obtain 
he ee Sih ae ke Ae (4.22) 
which leads directly to the equation 
pe 7M(k-N-1) = r"(k-N- 1} — Kysie™(k) (42) 


Equations (4.21) and (4.23) define the Lattice form of the whitening filter (see 
Chapter 3, Section D). This is also sometimes referred to as the analysis model. 


The structure is illustrated in Figure 4.1. 
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In order to develop a synthesis model we need only to solve equation 
(4.21) for e%(k) since we know that e%(k) is equivalent to y(k). Therefore, the 


synthesis equations are 


eN(k) = eNt(k) + Kyy yr (k-N-1) (4.24a) 


r+t(k—N—1) = rN(k—N—-1) — Kyyie*(k) (4.24b) 


The resulting structure is as illustrated in Figure 4.2. Note that in this model, 
y(k) is indeed being expressed as a weighted sum of the backwards errors, which 
represent an orthogonal coordinate system. 


This concludes the discussion of the classical Lattice formulation. 


B. GENERALIZED ORDER UPDATE RECURSIONS 
In this section a more general] linear prediction problem is considered. No 
assumptions are made as to the origin of the data. In fact, the data need not 
represent a time series at all, and certainly shift invariance is not required. The 
ordering of the data is simply chosen in some convenient fashion. The generality 
of this formulation will allow its application to multidimensional and nonlinear 
problems. 
1. Definitions and Formulation 
In this section quantities are defined that will be required to complete 
the statements and proofs contained in the remainder of the chapter. 
A realization of the random process Y is given by the vector 
Beeo< A < Kk. 
The error. e,\,. in predicting the element y**! from the previous N 


elements of Y is given by 


eos hy (ke iy? for A = 0..... K (4.25) 
where 
meee 0° 0 hp Shy x... -:° chifl 0 --- ol (4.26) 
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r® (kK) ri (k-1) re(k-e) 


Figure 4.1: 1-D Lattice Analysis Model. 


e% (kK) 


r (ae 





Figure 4.2: 1-D Lattice Synthesis Model. 
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The norm of this prediction error is given by 


1/2 
| leva} | = | e(etsrh| (4.27) 


A normalized version of the forward error is given by 


N 
G1 = -- aie 
| | x41 | | 
= oS 2ke (4.28b) 
where 
: 1 
ay\(k+1) = a hy (k~+1) (4.29) 
| k+1) | 


\ ao: , 
The backwards prediction error, ry, is the error associated with the 


prediction of y*"% given the next N elements of Y. It is given by, 


men = hy(k-N)y? (4.30) 
where 
eek N)i= 0 --- OC 1 =hy ns, —hpinee °°: —he 0 --- 0; (esa) 


Then 
Ty-N = 3 xn (4.32a) 
rN 
= byS(k—N)y? ee 2b) 
where 
| (kN 
by ale : (4.33) 
Ty | 
and 


om, 
Lee 
Fg 
| 
i 


ony | | (4.34) 
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In order to generalize the results of Section A of this chapter we need to 
introduce two different families of coordinate systems. We will refer to members 
of these families as either forward or backward local coordinates. The term local 
is used because a different coordinate system will be associated with every value 
of k and N. We define the (k+1) - indexed, N-th order forward coordinate 


system as 


ue 


yo (kN le (4533) 


K 
¥ 


The corresponding coordinate transformation from the unprimed coordinate 


system to this local forward coordinate svstem is given by 





ay? (k—1.N) | : 


Oy’ 
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| 1 ae hens - = ea he vs | 

| 0 1 —hyws "Gael hea ah a | 

lo 0 1 | ae es a | 

| | 7 (4.36) 
O : 0 

| | 

| 0 1 alii 

lo - Ee ee 
ae 0 bal 


where the 0’s are zero matrices and the I’s represent identity matrices. 
Similarly, we define the (k-N)-indexed, N-th order backward coordinate 


system as 


ik N.N) = | - (4.37] 


_k N-1..k-1 Na k= Neel 
Ve ihe y Se py 


The corresponding coordinate transformation is given by 


oe (k- N.N) 


dy? 
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Uo. oe ee 
| 1 0 0 0 | 
—hy wes 1 0 0 0 | 
| —hyine: —hyen+e J 0 0 | 
"I ale (4.38) 
| - | 
| —hXN., hee —hewis °° 1 oF 
JT aba chi ohvtNs 2 cheat a | 
Ol 0 I 


These definitions are sufficient to state and prove the theorems presented 

in the remainder of this chapter. 
2. Generalized Levinson Algorithm 

The Levinson algorithm (equation (4.14)) can be extended to recursively 
compute the forward and backward prediction coefficients defined in equations 
(4.26) and (4.31). In this section two forms of the algorithm are presented and 
proved. First a non-normalized version is introduced, then using this result a 
normalized algorithm is developed. 

a. Theorem 4.2: Generalized Levinson Recursion (Regular Form) 

The forward and backward prediction coefficients defined in 


equations (4.26) and (4.31) can be updated using the following recursive 


equations. 
oan : ey 
fhe ee keel oe ee D - (4.39a) 
hos i 
N | | 
Pa : one 
eee k= Nes Oe (4.39b) 
Chey, | 
where 


ree = efraants} (4.40) 


b. Proof 
The forward and backward prediction errors are scalars so their 


representation is identical in all frames of reference. Therefore. we can write 
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eX, = hN(k+1)y? (4.41a) 


= K§(k+1)y* (k+1,N) (4.41b) 
and 
FN = hN(k—N)y? (4.42a) 
= KN-(k-N)y*” (k-N,N) | (4.42b) 
where 
Peri = bY (ken)? — (4.43a) 
dy* (k+1,N) 
ae 7 O —Kyr yoy on -K, iP () os 0) (4.43b) 
. 
hN(k+1) = KN (k+1)22 aot (4.44) 
Oy 
and 
Ky-(k-N) = i’ (k-N}——2 (4.45a) 
dy* (k—N.N) 
ae Oeste, = —-K, 0 --- 0 (4.45b) 
. y EEEAS Sf tec. 
hy (k—N) = KN (k-N)o2 SAND (4.46) 
V 


The normal equations in the primed and unprimed coordinate 
systems can be solved for Kf\(k—1) and K,S-(k—N). This is once again (see Section 
A of this chapter) a straight forward matter because the correlation matrices. 


Ey? (k+1.N)y* (k+1,N)}] and [Ef{y? (k-N.N)y* (k-N.N)}]. will contain only 


diagonal elements. The solutions are 


Pek 1) = 
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efi uray} ebiay'-1ai] 


efor acsisyrh efota-u.xy} 


Giecere 150) >= =30 


(4.47) 
Ky\-(k-N) = 
Eftk-Nyt0-N] efy(k-NMA-NN 
O° 8+ 0 aa. 6 eee) 
Bfot any efoto-nnn 
(4.48) 
Consider the (k-N)-th term of the (N+1) order version of (4.47). 
Ebi uy aeaay} 
‘eS ——————— oe (4.49a) 
E{ot so.) 
efvik-nnts| 
——_--——— (4.49b) 
E{int wr} 
Efe in N i 
are ane (4.49c) 
clint 7} 
nee k 
N PN-=1 (4.49d) 
T-N 
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Similarly, 


or N a 
F(k—N) = Saree aN (4.50) 
eK ey | | 


For the N=1 model the forward error prediction coefficients are 


(from equation (4.44)) 


A 
Ihi(k+1)) = K}(k+1) 27 GALS) (4.512) 
y 
=f0--- 0 =K 10--- 0 (4.51b) 
O | 
= 0 oe () to 10-:°: 0) (4.51c) 
fom: | 


Equation (4.39a) yields 


a | | ey. 1 5 
hd(k+1)) = [h2(k+1)] — py [h?(k)}— ron (4.52a) 
| k 
| 0 
= {0 010 | — p#- ski , 0 010 0 (4.52b) 
[ecllei et 
’ 0 ; 
| e ' : 
mere 0 = py . ~~ | tas). ~ race (4.52c) 
esse 


Therefore. the recursion of equation (4.39a) holds for N=1. We assume it is valid 


for N and verify it for N+1. From equation (4.44) we have 


et —1) = KX“ (k+ 1} 


Pee 7 
(k= 1.N) e 
—_ (4.53a) 
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0 
—~ KAY (k+1) 
— Ke yea(k+1) + ona Kill (k+ 1) 


ale (4.53b) 
= KN ee hi Keke eR (ee 
1 
0 
0 
= ihN(k+1), = KP v(keO 0 1) heen (4.53c) 
at ae 
= hX(k-1) - ———— px. hi (k-N) (4.534) 
Th oN 


This completes the proof of equation (4.39a). Using similar 
arguments the backwards recursion of equation (4.39b) can be verified. 
c. Theorem 4.3: Generalized Levinson Recursion (Normalized Form) 
The normalized prediction error coefficients defined in equations 


(4.30) and (4.34) can be recursively updated using the following recurence 














relations: 
ee) \ 
a, (k—1) a;'(k=1) 
| Say 4.54 
fee L (5 ‘S O(pr_,) WAUSS ( 
where 
: | ] aN 1 a 
a | (4.55) 
V eee No 


and px., is given by equation (4.40). 
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d. Proof 


The proof follows directly from the non- normalized version of the 
generalized Levinson algorithm (equations (4.39)) and the definitions of the 
normalized prediction error coefficients (equations (4.29) and (4.33)). 


Using (4.29) in (4.39a) we obtain 


aN*1(k+1) 

1 | | ex. | | : 

Sal | deta | aX (e443) - pt SO | | ni] | biY(kK-N)] (4.56a) 
resi | | | | ren | ; 
Dalee ert | Ni}, k bk N ASSGh 

Memes 1) — prsaba (kN) (4.5 ) 
lt +1 
Similarly. using (4.33) in (4.39b) we obtain 
N+1 : itn | N k .N Sue 
pees N) = ————-—— [by (k—N) — pyiiay (k+1), (4.57) 
LiMn | | 
The proof of the fact that 
Mee) & end 1 (4.58) 


1» .N+1 to N41 ‘ 
ey 1 | Tk-N V1 = (px)? 

is relegated to the next section. 
The initial values for forward and backward normalized 
autoregressive parameters are obtained by setting N=O in equations (4.29) and 


(4.33). This yields 


ay(k+1) 0 ae 0 ee Cee ae a (4.59a) 


by{k) = 0°--- 0 —*"—— 0: :- 0 (4.59)) 


We note the similarities in the two generalized forms of the Levinson 
algorithm presented in this section with that presented earlier in this chapter. A 
little reflection will convince that equations (4.14) are simply a special case of 


equations (4.39). 
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3. Error Order Update Recursions 
a. Theorem 4.4 
The (N+1) order errors can be calculated from the N-th order errors 


through the recursion 


Cea Cet 
+ + 
=yir | = lens) | on (4.60) 
—N I KN 


where O(px,,) is defined in equation (4.55). 
b. Proof (Outline) 
Using the normalized form of the Levinson algorithm (equation 


(4.54)) and the following relationships 


Gri = aN (k+1)y? (4.28b) 


eo) = ay (ke L}y? 


Toe be ky (4.32b) 
Ton = by *'(k-N)y? 


equation (4.60) is easily verified. 


We now return to the proof of equation (4.58). That is. 


NN - 
Cxaim | Tew 1. 1 


7 ONS a = 5 = aS —_— 4.58 
ena | | Mon | </ 1 = pee Ge) 


Working with the term on the left-hand side of the equation 


‘#*" 


a ieee: 
E ef 
neice 7 


ie = hese — 
SVE +1 f 


ee (4.61b) 
yest 
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kt 1—N , N 
rhasBh ah ea (4.61c) 


Ve (pNa1)° 
(4.61d) 


Similar arguments can be used to verify the other relationship given in equation 
(4.58). 

Figure 4.3 symbolizes a single stage of the recursions of equation 
(4.60) while Figure 4.4 illustrates a thirl order analysis model. 

The interesting feature of equations (4.60) is that they do not make 
any assumptions about the nature of the given data. The data values need not 
be delayed versions of each other, as is the case for the autoregressive model. 
Any set of data values can be used. This fact will be significant when we deal 
with 2-D and nonlinear lattices. 

4. The Generalized Schur Algorithm 

In this section an algorithm will be presented which allows the 
calculation of the partial correlation coefficients in a direct manner. It will he 
shown that knowledge of the correlation matrix is sufficient to calculate all the 
reflection factors and thus solve the normal equations (by use of the Levinson 
algorithm). The method used has come to be known as the Schur algorithm Ref, 
an 

In order to obtain the desired result we must introduce two new 


quantities. defined as 


97 





(Cg) pee 






e ae 
j2 
—~Onar /{1-[Oned 3 
1/2 
+ One / ft=~[Oned 17 
_ ie 


a. Single Stage of Generalized Lattice Filter 
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b. Alternate Representation 


Figure 4.3: Generalized Lattice Filter Sections 
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Figure 4.4: 3-rd Order Generalized Lattice Filter 
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ax... (k+1) = eet. SP kes aria (4.62) 


Bx+i(k—N) = eft} = b,(k-N)R” (4.63) 


where 


R#* = ofr (4.64) 


is a correlation matrix. 
a. Theorem 4.5: Schur Recursions 
The quantities ax(k+1) and #x(k—N) defined in equations (4.62) and 


(4.63) can be updated according to 





A A 
ax.,(k+1) | ax(k+1) | 
. | = Ope, ; 4.69 
ai(k—N) (px) 3x(k-N) | 
and the partial correlation coefficient. px,,, can be calculated from 
k-N 
k On (k+1) 
16 4.66 
PN +1 Bat *(k—N) ( ) 
bor yoo 


The proof of equations (4.65) follows directly from equations (4.60) 
and the above definitions for ag(k=1) and 6.2(k—N). 
Beginning with the definition of the partial correlation coefficient 


given in (4.40). the relationship given in (4.66) is verified as follows 


xe efecanta} (4.67a) 
efedat}esiom (4.67) 


2 rfetat Stitt -N) (4.67¢) 
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= an N(k+1)byin(k—N) (4.67d) 


it 
1 1 

by n(k—N) = Ty) 1) > TRIERIENTT? (4.68) 
and 

Bem s—N) = [RUONKCN))}/ > (4.69) 
Therefore, 

Stel a 
mass (4.70) 


Be (kN) 


Using the initial conditions given by (4.59) the following initial 


values are determined for the Schur algorithm 


1 


ag'(k+1) = es Rea, (4.7 1a) 
ad(k-1)! = WE ssa ale Rie... Roe (4.71b) 
Pek) = a Re (4.72a) 
33(k) = a" el Ra... RAK (4.72b) 


where the parenthesis used to surround the first indices in equation (4.71) simply 
indicate that they are fixed at the value indicated. 
The Schur algorithm implies a filter structure identical to that of 
Figure 4.4. In this case the input vectors are the rows of the correlation matrix 
(normalized by the square root of the diagonal elements). 
5. Synthesis Model 
The original data. Y. can be svnthesised from the model parameters 


obtained from the Levinson algorithm. equation (4.54). It is also possible to 


regenerate the y* directly from the lattice parameters. The desired result is 
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obtained by solving equations (4.60) for &., and Avy. 


—N —N —N ou 
ex+1 = 1 — (pN41)” ear + PX-Tk_N (4.73a) 
<N —N —N+1 

Tk-N = 1 — (enai)? Thon — PN+1eKs1 (4.73b) 


Equations (4.73) constitute the synthesis model. They imply a structure 
similar to the analysis model of Figure 4.4. but with the direction of flow of the 
forward error signals reversed. The processing performed at each stage can be 
visualized as in Figure 4.5. A complete third order synthesis model is shown in 
Figure 4.6. Each horizontal path in Figure 4.6 represents a separate synthesis 
model. synthesising a different component of Y. The coordinate system for each 
of these models depends only on values of y* which appear farther down, that is 
thev have a smaller value of the index, A. 

Compare Figure 4.3b and 4.5b. It is apparent that the behaviour of the 
backwards error signals is identical in the two cases. Hence, it is possible to 
construct a synthesis model that only reverses the direction of the forward error 
corresponding to the point being predicted. This assumes knowledge of the other 
inputs (zero order forward errors) to the lattice. Such a configuration is 
illustrated in Figure 4.7. 

The amount of knowledge possessed about the signals used for the 
predictions dictates which form of the synthesis model should be used. If little is 
known. then estimates must first be generated which can then be used in the 
prediction. This corresponds to the model of Figure 4.6. If complete knowledge 
is available (either from initial conditions or previous predictions) then Figure 4.7 
can be used. It is also possible to construct models which exploit partial 
knowledge of the input signals and thus fall between these two extremes. In this 
case the known signals should be input as zero order forward errors while the 


unknown ones must be estimated. 
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6. Stochastic Fourier Series Interpretation 


From equations (4.72) and (4.40) we can deduce the following result 


N 
ee = eee a > efetaintalr (4.74a) 
d=1 
N+] : d ee 
= ea, + yy Pcie | Px eit a (4.74b) 
\=1 


where e., is equivalent to y**!. These expressions offer an alternate interpretation 
of the lattice filter. Equations (4.74) describe a stochastic Fourier series 
expansion of the forward error sequence where the basis functions are given by 
the backwards error signals. The Fourier coefficients are related to the partial 
correlation coefficients. 

This concludes the review of existing lattice formulations. In the next 
chapter new. multidimensional extensions of this theory are presented. In 


Chapter VI these results are used to derive original nonlinear lattice structures. 
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Figure 4.5: Generalized Lattice Synthesis Filter Sections 
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Figure 4.6: 3-rd Order Generalized Lattice Svnthesis Filter 
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Figure 4.7: 3-rd Order Generalized Lattice Synthesis Filter 
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Vv Ewe DIMENSIONAL LATTICE STRUCTURES 


In this chapter a new two-dimensional lattice structure will be derived and 
discussed. Lattice modelling of two dimensional fields has recently received 
considerable attention [Ref. 43,44,45]. In one dimensional lattice modelling, 
updating the order introduces only one new point into the model support. An 
order update in a 2-D lattice model must introduce O(N) new points, where N is 
the model order. Several different solutions have been suggested to this problem. 
The first due to Marzetta [Ref. 43,44], uses a particular ordering of the data to 
reduce the problem to one dimension. He proposed a half-plane support which is 
infinite in one of the two dimensions. This approach, while maintaining several 
of the nice characteristics of 1-D lattices. such as correlation matching and 
producing a minimum phase filter, leads to very long delay filters. 

A different approach, proposed by Parker and Kayran Ref. 43). 
simultaneously introduces many points into the support when the model order is 
increased. Their filter uses a quarter plane support and introduces three 
parameters at each order update. Therefore. it lacks sufficient parameters to 
represent all classes of 2-D autoregressive quarter plane filters. More 
importantly. it lacks the property of orthogonality so that the cascading of stages 
does not lead to an optimum filter (better filters are possible using an equivalent 
number of parameters). It’s simplicity is attractive and good results have been 
reported using this approach [Ref. 46]. 

The theory presented here maintains features of both previous approaches. It 
utilizes the generalized lattice theorv presented in Chapter 4 to decompose thi 
global O(N) point update into O(N*) single point local updates. It maintains the 
portant property of orthogonality so that the solution at all stages is optimum. 
Although only the quarter plane support case is presented here in detail. the 
theory can be used for any shaped support. It is shown that the Levinson and 


Schur algorithms (see Chapter 4) can be used to solve the 2-D linear prediction 
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problem. In its most general form the lattice contains O(N*) parameters while 
there are only O(N?) points in the filter support. Several structures are presented 


which take advantage of shift-invariance and reduce this requirement to O(N’). 


A. GENERAL FORM OF 2-D LATTICE FILTER 

The theory used is exactly that presented in the previous chapter. The 2-D 
structure results from a careful selection of input data. To illustrate the 
proposed 2-D lattice structure we will consider a 2-D linear prediction problem 


which utilizes a quarter plane support. The 2-D data field is given by 
Y= a = vrs (5.1) 


where apap.) cuelee=mlvam 


where L* is an index set given by 


de < x} (a4) 


Points will be used from this data field in a particular, convenient order. We 


define an ordered index set 


i {0.04 (1,0),(0.1420),02, ay (N,0),(0,N). 


(it) (2en)(3.2). 2, (51) eee 


\ 


(K—2,K),(K—1,K- 1).(K.K-1),(K- LK) KK) (5.3) 


Other orderings are possible and equally valid. This one is chosen merely to 
illustrate the concepts. The desired ordering of *“L* . to obtain 6L® . can be 


accomplished by the following. computationally efficient algorithm 
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k= 0 
ror m= 0 to 
for nO = 0 to 2*m0 
if (mod(n0.2) = 0) 


then begin 
i= m0 
j=n0/2 
end 
else begin 
i= n0/2 
j) = m0 
end 
BLK (k) = *LX (iJ) 
k=k+1 
next nO 
next m0 


In this algorithm LX (k) has been used to describe the k-th element of the index 
set éL*. while *L*(i,j) has been used to indicate the (ij)-th element of *L*. The 
(K+1)? elements of 7L* have been ordered into a one dimensional index set. éL*. 
The elements of éL* can be numbered. consecutively, from 0 to (K+1)*-1. The 
notation (k,l) - q will be used to mean: the element of the index set corresponding 
to the q-th element prior to the element (k.]). For example, (2,0)-3 would 
indicate the element (1,0) (see equation (5.3)). Occasionally this notation will be 
abbreviated to simply, kl-q. 

Define the (q-1) order. normalized. forward error associated with the 
prediction of the point y(k.l) from the previous (in the sense of §L*) (q-1) points. 


as 


a. i xe ee 
= avy (k.Dy *? (3.4) 


where the implied summation over (A,.A,) < *L*. can be carried out in any order. 
as long as all components are considered. It is preferred to think of (A,.A.) 
K 


belonging to *L* rather than ¢L* as this maintains the two-dimensional character 


of the problem. 


The a}y\(k.l) can be interpreted as the components of a second order 


covariant tensor. These components. for a range of indices (A,.A2) < (k.J)--q (in the 
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sense of (5.3)), or for indices (A,.A2) > (k.l). are equal to zero. When (A,,A3) = (k,l). 


the component 


a (les = ama (5.5) 


where e797! is an unnormalized version of @3’. 


A normalized backward error associated with the prediction of y((k,l)-q) from 
the next q-1 points of gL*, is defined by 


—_g= on A,A 
Fog = bAAU((K.)—q)y 1? (5.6) 


As with the forward error prediction coefficients, the byj)((k,l)-q) can be 
interpreted as components of a second order covariant tensor. The components 
b#y((k-l)-q) equal zero for the range of indices (A,,A2) < ((k,1)—q) or (Ay.A2) > (k,}). 
For the case when the index (A,,A.) = ((k,l)-q) the component 


_ I 
bf ((k.1)—q) = = 
a a | 
where r,j_} is an unnormalized version of r774. 
1. Normalized 2-D Levinson Algorithm 
a. Theorem 5.1 


The 2-D prediction error coefficients can be updated in order using 


the following recursions 


ayy (kl) ata (k.l) 3) 
= © k} | on 
be. (Uk) a) | > 8? | gste(k.t) a) 
where 
’ I ae ; 
OlP, Vian KI)? ie J (5 °) 
and 
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b. Proof (Outline) 

The proof follows directly if the index (k.l) is replaced by a single 
index which runs from O to (K+1)*--1 and thus indexes the elements of 6L*. The 
equations (5.8) are identical in form to equations (4.54) and the proof presented 
there can be applied. This approach is equivalent to reordering the data field 
into a vector in the order specified by 6L* 

An alternate proof is possible by peerolinine the methods of Chapter 
4. Alternate coordinate systems could be introduced (similar to (4.35) and 
(4.37)). The necessary transformations can be found and all the steps of the 
proof of Theorem 4.2 can be generalized to these higher order objects. These 
concepts will not be explored further here except to note that they could be 
extented to solve problems in any number of dimensions. they are not restricted 
to the two dimensional case being studied here. 

2. Normalized 2-D Error Order Updates 
a. Theorem 5.2 
The two-dimensional prediction errors can be updated in order 


according to the following equation 


ao asc 

&x) "5 k] - 

—q =o) low) 
Tykl-g kl—q 


b. Proof (Outline) 

The proof follows from equations (5.8). If an inner product is formed 
on both sides of the equation (5.8a) with the data field Y equation (5.11a) 
results. Similar arguments can be employed in the verification of (5.11b). 

2 


ee 2-)) Form of Schur Recursion 


Define the quantities 


ra,\ Laghyd : 
oh es) alae (5.12a) 


An 
my (k.)) = eles 'y 


8,1 (k,l)—q) = en” nt = bt NE oye aa (5.12b) 
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where 


eae = efi (Serer 


We note that for the 2-D case the correlation is a fourth order tensor. 

The Schur recursions for the 2-D case are then given by the following 
theorem. 

a. Theorem 5.3 


The generalized 2-D Schur recursions are given by 


A,A A,X 
a, * (k,l) 0 6 27?(k,1) 


= O(p x) (5.14) 
By *((k.l)-4) | Beat 7C(1) 4) 
and p,' can be calculated from 
eke 
pi = @ gai (ki) (Salo) 


By ((k.1)—q) 
b. Proof (Outline) 
The proof follows identical arguments to those of Theorem 4.5 and so 
is not given here. 
4. 2-D Lattice Structures 
The aenvation® presented do not assume shift- invariance. Models are 
built for each point in the data field starting with the point y(K.K) and ending 
with the point y(0.0). All the models are not equivalent. however. In fact. no 
two models are identical. The only model that is quarter plane is the first one. 
that is the model corresponding to the point v(K.K). A quarter plane model for 
any point. y(m.n), in the field can be built by considering an appropriate subset 
of the set Y (equation (5.1)). The subset would start with the point y(m-N.n-N) 
and continue in a quarter plane manner until the point y({m.n). for an N-th 


(elobal) order immodel. This support can be written 


j= fal (5.16) 


where 


eZ 


where 

i = {im-N)m-=1) oe a} (5.17a) 
and 

i = {-Nhia-Ne), a ; (5.17b) 


The (N+1)? terms of this index set can be ordered to form 


aN = {im-Nan-S}.m-N- 1,n—N),(m—N,n—N-+1)..... 


(m-tan-t}mnn- 1} (m1). (5.18) 


The subset of Y given by (5.16) and the ordering of (5.18) are illustrated in 
Figure 5.1. This could be done for each point in the entire data field Y. yielding 
(K+1)? models. If the process is shift-invariant. then all the models would be 
identical. Other simplifications in the model are possible if shift-invariance can 
be assumed. These will be the topic of the next section. In addition. if 
ergodicity is assumed the required statistical averages can be replaced by 
appropriate spatial ones. 

A second order quarter plane 2-D lattice model] for generating the 
prediction error associated with an arbitrary point y(m.n) is illustrated in Figure 
5.2. In this diagram the forward and backward error fields are indicated 
pictorially rather than symbolically. The icons used are defined in the legend on 
the diagram. The large squares, at each stage. indicate the entire support for the 
second order model. The small blank squares indicate the additional support 
(besides the error field squares) used to generate the given error fields. For 
example, the forward error field in the upper right hand square is generated by 
predicting y(m.n) from all the remaining data points in the large square. 


including the one indicated as a backward prediction error. 
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Figure 5.1: Filter Support Ordering (see equation 5.18). 
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Ds 


] 


Figure 5.2: Genera 


The doubly hatched squares. corresponding in this diagram to zero order 
errors, are inputs to the filter. In general, (in later diagrams) doubly hatched 
squares are considered to be inputs. although. they may not always correspond to 
zero order errors. 

The ordering chosen for 6L* (equation (5.3)) is only one of many that 
could have been chosen to implement a quarter plane model. This particular 
selection was made for ease of implementation and because it guaranteed that for 
every 2N+1 local updates. a global order update would be completed. In Figure 
5.2 the filter can be visualized as a cascade of increasing order filter sections. For 
every global update. O(N*) local order updates must be completed. This implies a 
total of O(N‘) updates for an N-th order filter. In general this is too large a 
number to allow these filters to be used for any real time applications. 

In the next section the problem of reducing the complexity of the 2-D 


lattice filter is examined and some solutions are proposed. 


B. REDUCED COMPLEX TTY 2-Dar a Gea iis 

The assumption of shift-invariance allows certain of the backward prediction 
errors to be considered to be shifted versions of each other. This eliminates some 
calculations. Various structures are possible depending on the types of shifts 
introduced. We note that no single type of shift (neither 2, !.2.),2,!'2."') will 
introduce all the new data elements into the support that are required for a 
global order update. Because of this. additional prediction errors will have to be 
introduced at each stage. This reduces somewhat the advantage gained by the 
shift-invariance assumption. 

Two types of delays will be considered in detail. Initially. several models 
involving diagonal] shifts are examined. Later. a model involving a horizontal 
shift is discussed. We begin by introducing a diagonal shift operator. D. This is 


=| 


equivalent to multiplication by z, ‘2. ' in the bivariate z-transform domain. 


~ 


Because of the assumed shift invarianee the following statements can be 


made 
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R(m.n,m-—i.n—j) 


ef(m.aiyim-in-ah (5.19a) 


= p{Div(ma)Dsim-in-))} (5.19b} 
= y(m—tn-t)y(m-i-ta-j-n)} (5.19c) 
= R(i,j) (5.19d) 


where R() is a correlation function. The correlation is only a function of the 
relative positions of the two points not their absolute positions. If we adopt for a 
moment the Hilbert space formulation (see Appendix A) we conclude that the 
diagonal shift operator is an inner-product preserving operator and so the use of 
the shifted versions of the backward error signals is justified. 

Consider the structure illustrated in Figure 5.3. It represents a third order 
quarter plane lattice filter. At each global stage, (2N-1)°-3. lattice coefficients 
are introduced. Therefore. an N-th order model requires O(N*) parameters. 
Notice that at each stage two new errors are introduced. They each require the 
solution of an (N-1)° point prediction problem. For small N this is an 
insignificant number, however, for large N it becomes overwhelming and the 
required number of parameters again becomes O(N‘). It is difficult to analyze this 
structure analytically. as the index sets for each prediction error are different. 
The support for different errors follow different patterns. This. and the 
complexity issue make it a structure that is really only of academic interest. 

The structure of Figure 5.4 is a true O(N*) parameter inodel. It avoids the 
addition of the new error signals at each stage by introducing them at the outset. 
The structure has a support that differs slehtly from quarter plane. A more 
significant drawback. however. is that the inaxinum order of the filter must be 
fixed at the start. If the maximum order is overestimated then some unnecessary 


computations will have been performed. If on the other hand. the maximum 
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Figure 5.3: Reduced Complexity, 2-D Quarter Plane Lattice Filter. 
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required order was underestimated, then a great price must be paid to increase it. 
However, this is considered to be a superior structure because of the regularity 
and complexity reduction it offers. Therefore, a more detailed analysis of this 
model will be performed. 

The ordered index set (equivalent to (5.18)) in this case (for N=2) is given 
by 


oe. = {m-2n-2) (m=1.0-8)(m—8n—1}(m=19=2),(- 29-2) 


(mtn 2}(mBhm—2a) (sn—1},m-1.) (5.20) 


The following relations can then be used to advantage to update the backwards 


errors 
ey e-1) = Dita) ou) 
ee 2-1) = Ditin-1n) (5.21b) 
<a oe Din aa ae) 
7 ae bane (5.21d) 
ee = Dia a (5.21e) 
eee n-2) = DT(ma0-1) ee) 


In general. whenever (m-i.n-j) equals (m.n)-q. then 
et ss 0) = Om -i 3) (5.22) 


Equation (5.22) is a simple rule for exploiting the shift- invariance of the data 
field. 

One final reduced complexity lattice model will be introduced. It will serve 
to illustrate the variety of structures possible and in particular will yield a model 
for which an especially convenient synthesis model can be constructed. The 


model incorporates a horizontal (z.') shift. rather than the diagonal shift used in 
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Ficure 5.4: Reduced Complexity. O(N‘), 2-D Quarter Plane Lattice Filter 


the previous models. A horizontal shift also is a inner product preserving 
Mapping so that its use 1s permitted. 
The ordered index set used for this model is given by (for the second order 


case) 


Abe. = {im-2.n-2)(m-1 -2),(m,n—2),(m—2,n-1), 


fm-tan-1},{mn-1)(m-2a},m-1.8}, (ma) (5.23) 


Define a horizontal shift operator, H, by the relation 
y(m.n—1) = Hly(m,n)| (5.24) 


The following relations hold due the assumed shift- invariance of the given data 


field 


fee) = Hirinn)! (5.25a) 
ee 1) = H Fn.) (5.25b) 
= 2.n- 1) = 151 Fee (5.25c) 
Pee 1) = Etisal (5.25d) 
Beans = Hina), | (5.25e) 
Wee 1) = Hotes.) fomeor) 


In general. whenever, (m-i,n-j) equals (m,n)-q then 


= 4 q 
Mim=i.n- 53-2) H Pim -a.n- 


nt [oa6 | 


Using these simplifving relations and equations (5.11). the model of Figure 
5.9 can be deduced. This model still contains O(N*) parameters. however. the 


actual number is only about one quarter of that required by the previous 


structure (Figure 5.4.) 


PAM | 


This algorithm shares with the previous algorithm the shortcoming of having 
to estimate the maximum order of the filter at the outset. Despite this. it is 
believed that this model offers a good compromise (probably the best to date) 
between model accuracy and implementation complexity. In the next section it is 
demonstrated that the synthesis form of this algorithm has some particularly 
desirable properties. In Chapter 7, it will be shown that this algorithm is well 
suited for highly parallel VLSI implementation. 


C. SYNTHESIS MODEL 
The synthesis results of the previous chapter are easily extended to the 2-D 

case being considered. The data fields can be regenerated from a knowledge of 

the lattice coefficients and the forward error fields, it is not necessary to explicitly 

calculate the forward and backward error prediction coefficients. 

1 


The desired result is obtained by solving equation (5.11) for 43 


This yields 


att = /1— ob at + oti (52a 


and Gees 


Ti = i‘ [pave igigee jae ps ed (5.27b) 


These equations establish the method for regenerating the original data field. 
They describe the processing that must be carried out at each stage of the 
synthesis process. As an example of their application. consider the second order 
synthesis model pictured in Figure 5.6. It is the synthesis counterpart of the 
reduced complexity analysis model of Figure 5.5. In order to regenerate the 
original data field processing should be carried out horizontally by rows. For this 
second order model, two rows and two columns of initial conditions must be 
specified. The required zero order forward error sequences will always be 
available from either the given initial conditions or from previous estimates. 

It will be shown that for VLSI implementations it will be more convenient to 
estimate all the zero order forward error sequences. This necessitates that three 
residuals to be input and that all forward error channels be reversed. Such a 


structure is illustrated in Figure 5.7. 
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Figure 5.6: 2-D Quarter Plane Lattice Synthesis Model 
Involving Only a Single Input Residual 
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A synthesis model corresponding to the analysis model of Figure 5.4 would 
have to process the data along diagonals. The required zero order forward error 
sequences would not be available at each stage in this case and so only a model 
analogous to the one of Figure 5.7 can be used. This would require the provision 


of five input residual signals. 


D. SYSTOLIC IMPLEMENTATIONS 

In the past, most signal processing algorithms were implemented largely in 
software due to their high complexity. More recently, with the advent of VLSI 
technology, there has been a shift towards specialized hardware implementations. 
These offer increased performance at a low per-unit cost. A particularly 
promising class of implementations, first suggested by Kung and Leiserson [Ref. 
49]. are the so called systolic arrays. These attempt to partition the required 
computations in time and space over an array of identical processing elements. in 
order to increase throughput. 

Kung |Ref. 50: pp. 869] defines a systolic array as a network possessing the 


following features: 


a)Svnchrony: The data are rhythmically computed (timed by a global clock) 
and passed through the network. 


b)Regularity ie Modularity and Local Interconnections): The array should 
consist of modular processing units with regular and Sau e) local intercon- 


nections. Moreover. the computing network may be extended indefinitely. 


c)Tempora] Locality: There will be at least one unit-time delay allotted so that 
signal transactions from one node to the next can be completed. 


| Pipelinability (i.e. O(M) Execution-Time Speed-Up): A good measure for the 
efficiency of the array is the following 


Processing Time in a Single Processor 
SpeedUp Factor = ———S— 
Processing Time in the Array Processor 


A systolic array should exhibit. a linear-rate pipelinability. 1.e.. it should achieve 

an OU speceaveyy terms of processing rates. Where Nias the number of pro- 

cessor elements (PF‘s)}. 

Methods have been proposed for transformime algorithms into Systolic 
Mnplementations beginning with either an algorithinic description {sec 
Moldovan [Ref. 51]) or a signal-flow-graph (SFG) description (see Kung [Ref. 


50]). In this chapter we shall be using the second of these methods to transform 
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one of the 2-D lattice structures into systolic form. Rather than present a 
detailed discussion of the rules used in this method, we shall simply state them 
and then apply them to produce the desired systolic array. It is hoped that this 
example will clarify the procedures used. If further insight is desired the reader is 
referred to the cited reference. 
1. SFG Transformation Procedure 
The procedure used is based on a cut-set approach. According to Kung 


(Ref. 50, pp. 870] a cut-set is defined as: 


A cut-set in an SFG is a minimal set of edges which partitions the SFG into two 
parts. 


He proposes and proves that the following two rules that can be used to 


transform any computable SFG into a systolic array: 


Rule (i) Time-Scaling: All delays D may be scaled, i.e., D > D’, by a single po- 
sitive integer . Correspondingly, the input and output rates also have to be 
scaled by a factor (with respect to the new time unit D’) .... 


Rule (ii) Delay-Transfer: Given any cut-set of the SFG, we can group the edges 
of the cut-set. into "inbound edges" and "outbound edges", depending upon the 
directions assigned to the edges. Rule (ii) allows advancing k (D’) time units on 
all the outbound edges and delaying k time units on the inbound edges. and vice 
versa. It is clear that. for a (time-invariant) SFG, the general system behaviour 
is not affected because the effects of the lags and advances cancel each other in 
the overall timing. Note that the input-input and input-output timing relation- 
ships will also remain exactly the same only if they are located on the same side. 
Re icc they should be adjusted by a lag of —k time units or an advance of -k 
time units. 


2. Systolic Implementation of 2-D Lattice Filter 

Using the two rules given in the previous section the 2-D lattice synthesis 
mode! of Figure 5.7 will be mapped into a systolic array. There is some 
flexibility in the design. the result not being unique. The first choice that must 
be nade is that of the operation that is to be performed by the basic PE. In the 
case being considered several convenient choices are possible. The simplest 
element would be a multipher-adder. A slightly higher level operation would be 
that of the single lattice section given by Figure 4.5. A still higher level operation 
is conceivable hy grouping several of the lattice sections together. We shall use 
the second of these choices as it illustrates the general procedure without the 


added complexity inherent in the lower level implementation. 
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The mapping will be done in stages. Initially. the diagram of Figure 5.7 
will be redrawn in an SFG format. Rule (i) will be used to scale all the delays 
appropriately. Then, by successive application of Rule (ii). the SFG will be 
temporally localized (it is already spatially local.) The steps used are outlined as 


follows; 


(1) In Figure 5.8 the algorithm of Figure 5.7 has been redrawn in a SFG form. 
In this and subsequent diagrams delays are indicated by the letter D. The 
lattice section of Figure 4.5, as before, is indicated by the circles at the 


nodes. 


(2) Using rule(i), all delays are scaled by a factor of 6. That is, D ~ 6D’. This is 
indicated in Figure 5.9. The input and output signal rates must also be 


scaled by this factor of 6. 


(3) Using the cut-sets indicated in Figure 5.10 the delays are redistributed so 
that temporal locality is achieved. The resulting SFG is indicated in Figure 
5.11. Until now the processing at each node was assumed not to take any 
time. The delays going into each node can be combined with the lattice 
sections and be used to account for the processing time. In this way. the 
structure will appear as in Figure 5.12. In this last figure. the nodes have 
been shaded to indicate that the operation being performed within them 
consumes one time unit. 

3. Additional Remarks 
In this section we have shown that the 2-D Lattice structures derived in 
this chapter are ainenable to a systolic implementation. This is significant as the 
processing of 2-D data fields such as images in real-time requires high data rates. 

These rates can only be achieved in practice through the use of super-computers 

or specialized hardware. Due to the high cost of super-computers the second 

alternative is tle more practical. With the costs of VLSI production rapidly 
decreasing. it 1s now cost effective to produce dedicated chips even in very small 
quantities. For large scale productions the cost can be amortized over a large 


number of chips. vielding a low per-unit cost. 
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Although only a single specific implementation has been presented here, 
it indicates the ease with which the other algorithms discussed in this thesis may 


be transformed into forms which can be efficiently implemented in silicon. 


E. SIMULATION RESULTS 
The theory has been proven by computer simulation. Two different order 
models were excited by a white noise process. Several different order estimates of 
each spectrum were generated and compared to the originals. 
1. Example 5.1 
The first model simulated was described by 


.295y(m—1,n) — .470y(m,n—-1) + 0.0y(m—1,n—1) 


y(m,n) 


{ 


.055y(m—2,n) + .007y(m,n—2) + .003y(m—2,n- 1) 


— .O15y(m—1,n—2) + .022y(m—2,n—2) + u(m,n) (5.28) 


where u(m,n) was a 2-D zero-mean white noise process. 

The original spectrum is illustrated in Figure 5.13 while the first. second, 
third, and fourth order estimates are given in Figure 5.14. Notice that the 
original model is only second order so that the third and fourth order estimates 
can be used to examine the effects of over modelling. As can be seen from the 
figures the estimated spectra correspond very closely to the originals. Over 
modelling did not noticeably degrade the accuracy of the estimates. 

The actual algorithm used was that of Figure 5.2. Although it 1s 
unnecessarily complex. it is the most straight forward to implement. The 
generalized 2-D Schur algorithm was used to generate all the required lattice 
parameters. The computer subroutines used to accomplish this simulation are 
included in Appendix B. 

fe Example 3.2 
The second simulation used the following higher order autoregressive 


equation to generate the data 


y(m,n) = —.47y(m—1.n) — .03y(m.n-1}) + .195y(m-—1.n—1) 
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Figure 5.8: 2-D Lattice Filter of Figure 5.7 in SFG Form. 
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Figure 3.9: The Delays are Scaled by a Factor of Six. 
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Figure 5.10: SFG for 2-D Lattice Filter Showing Cut-Sets to be Used 
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Figure 5.11: Temporally Local Version of 2-D Lattice Filter 
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Figure 5.12: Temporally Local 2-D Lattice Filter. 
Shaded Nodes Indicate the Inclusion of a Delay 
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Figure 5.13: Original Spectrum For Example 3.1. 
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a. First Order Lattice Approximation of Spectrum of Examo!te 3.1. 
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Second Order Lattice Approximation of Spectrum of Example 5.1. 


Figure 5.14: Lattice Approximations Of Spectrum of Example 5.1 
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c. Third Order Lattice Approximation of Spectrum of Example 5.1. 
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d. Fourth Order Lattice Approximation of Spectrum of Example 3.1. 


.015y(m,n—2) + .055y(m—1,n—2) ~ .003y(m--2,.n—2) 


}- 


.0067y(m—38,n) — .015y(m,n—3) + .022y(m—2.n—3) 


ate 


.033y(m,n—4) — .085y(m—1,n—4) — .002y(m—4,n—2) 


— .000ly(m—2,n—4) + .0001y(m—4,n—4) + u(m,n) (5.29) 


where u(m,n) is a zero-mean white noise process. 

The spectrum of this model is shown in Figure 5.15. The first through 
fourth order estimates of the spectrum are shown in Figures 5.16. The general 
shape of the spectrum is identified in the first order model, although the fine 
detail is not introduced until the fourth order model. The position and relative 
magnitude of the peaks in the spectra are identified with great accuracy in the 
fourth order estimate. | 

In the next chapter lattice models are applied to the solution of the 


autoregressive nonlinear modelling problem. 
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Figure 5.15: Original Spectrum For Example 3.2. 


W 


SN fi 


re, eH | 
{/ _— * 
Pag lien” 7 






00 20 40 69 80 100 20 “40 1609 
HY ————————S— 


OM ALT ALD 
OSL LILA? 
KEILLOR LI 
Fea. os Og CBs TP a gas ose, thi 
Pn n= aetccscraccsstenceocesorereee Zl 
ses Perssszsz S25 Sera Og 
‘ae feel ay Spiga f 
ins sirsisrasessssstosseose OLY 
0 PIF t roe LOLOL L/D 
8 oy Qo scececzcosteceaceteete. 
5 ee -, ane 
0% ESS” 





a. First Order Lattice Approximation of Spectrum oi Example 3.2. 
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b. Second Order Lattice Approximation of Spectrum of Example 3.2. 


Figure 5.16: Lattice Approximations Of Spectrum of Example 3.2 
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c. Third Order Lattice Approximation of Spectrum of Example 35.2. 
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d. Fourth Order Lattice Approximation of Spectrum of Example 3.2. 
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VI. NONLINEAR LATTICE SIRUC TREES 


In the previous two chapters lattice structures have been examined in great 
length and have been applied to the problem of 2-D autoregressive modelling. In 
this chapter we return to the problem of nonlinear autoregressive modelling with 
the hope of solving the problem using a lattice structure. 

There have been several attempts in the literature to fit lattice filters to a 
Volterra like model. The first was due to Parker and Perry [Ref. 15]. They 
proposed a multichannel lattice implementation of an autoregressive-moving 
average nonlinear model. Their model was capable of providing an update in time 
of some terms of the expansion. It. however, does not introduce all the terms that 
are needed to constitute the next higher order model. 

A recent proposal by Zarzycki and Dewilde [Ref. 17] is a true generalization 
of the Levinson algorithm to the case of the autoregressive Volterra type model. 
In their work they considered a non-stationary nonlinear structure and showed 
that the stationary and linear models can be treated as special cases. Their 
model provides a true update in time order, introducing all the required terms. 
They found. as we did with the 2-D model. that not all terms could be 
recursively generated and that some would have to be introduced at each stage 
by other means. For these non-linear models only one error signal must be 
injected at each stage not two as in the 2-D case (if a triangular kernel is used.) 

The model proposed in this chapter is based on the alternate tensor form of 
the nonlinear model developed in Chapter 3. Recall that the model was based on 
interchanging the roles plaved by time order and nonlinear order. The lattice 
model presented here thus becomes recursive in nonlinear order. This means. for 
example. that the optimum cubic inode] is calculated from a knowledge of the 
optimum quadratic model. The theory is based on the generalized lattice 


concepts of Chapter 4 and the 2-D lattice structures developed in Chapter 5. 
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Although only the nonlinear order update model is discussed. some reflection will 
show that a time update can also be performed using the proposed model. 

For simplicity in illustrating the concepts only models involving two delays 
will be considered in this chapter. The theory is equally valid and easily extended 


to more complex situations. 


A. GENERAL NONLINEAR LATTICE MODEL 

As with the 2-D model of the previous chapter, the theory used is that of the 
generalized error order updates given in Chapter 4. A careful choice of the input 
data will yield the desired results. We begin by summarizing briefly the 
nonlinear model to be used, for the case of two delays. 


The estimate of the model output is given by (3.53) 


A New : 
mis= y (k-1) --- ye{k-N)H,, ..4, (6.1) 
To avoid unnecessarily complex algebra we consider the case when N=2. 
Aj Ay 
y(k) = y “(k-1)y “(k-2)Hy a, (6.2) 


The tensor product y(k-Ly"2(k—2) defines a second order tensor, Y(k), with 


components given by 


¥(k} = fy (k-Ny 2(k-2)] = 


y(k—2) yo (ke 2) = y)(k—2) 
y(k- 1) y(k-J)y(k-2) — y(k-A)y?(k-2) ++ -  y(k—3}y}(k—-2) 
BU) yk vy{k-2) y@ik—jyl(k—2) => yk—1)y®)(k—2) 
yiPi{k 1) ylPI(k is (k 2) vik —1)y lk a Bare y!P){k Hy") (k 2) 


(6.3) 


This tensor can be considered to be a shift-varying data field and can thus be 
modelled using the 2-D lattice model of Figure 5.2. A tensor of the form given in 


(6.3) can be formed for each time, k. If the process is time-varying then each of 
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these tensors must be separately modelled. If time-invariance can be assumed. 
then only one model need be developed. In that case the required correlations 
can be calculated over the ensemble of such tensors. 

A shift along the diagonal (or horizontal or vertical ) of the Y(k) tensor is not 
a inner-product preserving operation. This prevents the application of the 
reduced complexity algorithms of Chapter 5 to this nonlinear problem. 

The 2-D lattice model of the tensor Y(k) does not offer a complete solution. 
Notice that y(k), the sample that is to be estimated, does not appear in equation 
(6.3). Two possible solutions to this problem exist. First, the ordered index set. 
(5.3), can be extended by one element so that y(k) is included in the model. This 
has the effect of adding a channel to the general 2-D lattice structure of Figure 
5.2. A conceptually different solution is to first model the tensor Y(k) using the 
results of the previous chapter. The backwards error signals that result can be 
used as a basis for a Fourier expansion of the the sample y(k) (see equation 
(4.73) for details.) Both of these approaches lead to identical models but do 
provide alternate interpretations. The second approach will be the one used in 
the derivations of this chapter as it allows the results of the previous chapter to 
be apphed with little modification. 

The tensor, Y{k). can be considered to be a two- dimensional data field given 


by 


\ 


Y(k) = (a0 24 (6.4) 


Tenens = Le saetee 


where L? is an index set given by 


Ufo.) (6.5) 


Points will be used from this data field in a particular. convenient order. We 


define an ordered index set 
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6L? = {r..(0-t0)40-2) ---, (0,p),(p,9), 
(p—1,p—1),(p—2,p—1),(p—1,p—2), ee) (O.p—1),(p—1,9), 


-(41),0.0,(40),(00)} : (6.6) 


The (p+i)? elements of 7LP have been ordered into a one dimensional index 
set, GL’. The elements of gL? can be numbered, consecutively. from 0 to (p~+1)?-1. 

As in Chapter 5. the notation (m,n) - q will be used to mean; the element of 
the index set corresponding to the q-th element prior to the element (m,n). For 
example, (1.0) - 2 would indicate the element (1,1) (see equation (6.6)). This 
notation will often be abbreviated to simply mn-q. 

Define the (q-1) order, normalized, forward error associated with the 
prediction of the element y™(k—1)y"(k—2) from the previous (in the sense of éL’) 


(q-1) points, as 


aao) = af'(m.n)y (k-1y (k—2) (6.7) 


] 


where the implied summation over: (4,,A.) ¢ *L?, can be carried out in any order. 
as long as all components are considered. It is preferred to think of (A,.A,) 
belonging to “L’ rather than 6L? as this maintains the multi-dimensional 
character of the problem. 

The a 3(m.n) can be interpreted as the components of a second order 
covariant tensor. These components. for a range of indices (A,.A2) < mn—q (in the 
sense of (6.6)). or for indices (A,.A.) > (m.n}) are equal to zero. For the case when 


(A,.Ae} = (m.n). the component 


| 
<a (6.8) 


Ti fs 


ae (m.n) = 


where e,3,’ is an unnormalized version of e{.'. 
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A normalized backward error associated with the prediction of y((m,n)-q) 


from the next q-1 points of 6L*, is defined by 
Fgotg = bsai(mn—g)y (k—1)y %(k-2 6.9 
mas ¢) AjAo mn q)y ( yy ( ) ( . ) 


As with the forward error prediction coefficients, the b{\'(mn-q) can be 
interpreted as components of a second order covariant tensor. The components 
b#y,(mn—q) equal zero for the range of indices (A,,A2) < (mn—q) or (A,,A2) > (m,n). 
For the case when the index (\,,\,.) = (mn—q) the component 


b3,!,((m,n)—q) = ———— (6.10) 


| | tmtn—g! | 


where r@,/, is an unnormalized version of r2,",. 
1. Normalized Order Update Recursions 
The following two theorems are stated without proof. The proofs are 
identical to those of Theorems 5.1 and 5.2 of the previous chapter and so are not 
repeated. 
a. Theorem 6.1: Normalized Nonlinear Levinson Algorithm 
The nonlinear prediction error coefficients can be updated in order 


using the following recursions 


ayia,(m.n) a)’ (m,n) 
ue =a) | = Clectalaiaes (Guss) 
A y,{(m.n) q) x a,(™m.n) 
where 
] ~p mn 
On| =a : | fn | (6.12) 
vale aes a) Vas 
and 
Pics - frit | (6.13) 
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b. Theorem 6.2: Normalized Error Order Update Algorithm 
The nonlinear prediction errors can be updated in order according to 


the following equation 


Sinn Tao ' 
_ =Q(p.™) | __ (6.14 
a) et | at, 

In order to introduce the sample y(k) into the model we recognize 
the backwards errors as an alternate coordinate system. In particular the 


following vector is defined. 


_0 
T(0,0) 
T(0,0)-1 


*%) 


T(0,0)—2 


y= [Foo-a] = | (6.15) 


lS 
F(0,0)—(p+1)?+1 


Equations (6.14) correspond to a model structure identical with that 
given in Figure 6.2. In this diagram, the backwards errors given in equations 
(6.15) correspond to those that are shown leaving the uppermost stages of the 
filter. 


The forward error in predicting y(k) from the Y(k) tensor is given by 
Se = y(KieeKy-y* (6.16) 


Because the backward errors are uncorrelated it is a straight forward 


matter to find the (Fourier) coefficients. K,... They are given by 


hy < efi} Pppikln boa} _ efrbn ga ip} 


The m-th component is 
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o 
3 
| 


: eb Wreo-.} (6.180) 


Eero} (6.18b) 


\ 


| tex” | Pfs} (6.18c) 


[oe pelvars (6.18d) 
A normalized version of the forward error defined in (6.16) is thus given by 

te 1 zi ae ae 

am = —— ih) - SS obec ditier (ald (6.19) 


Recognizing that y(k) is a zero order forward error we can write a 


normalized version of the q-th order forward error as 


g—! 
= mo" waned eee Ie 
ee ee eee (6.20) 
k 


This last result follows if (6.19) is iterated. calculating successively 
higher order errors or from equation (4.73) where the equivalence of the lattice 
model and Fourier series was established. 

We conclude that expressing y(k) as a stochastic Fourier series (with 
the backward errors as a basis) amounts to nothing more than the addition of a 
supplementary channel to the existing lattice structure. Equations (6.11) and 
(6.14) hold for this additional channel and can thus be used to update the 
normalized error signals and inodel parameters. Figure 6.1 illustrates a quadratic 
nonlinear lattice filter. 

2. Uniqueness Of Lattice Parameters 
a. Theorem 6.3 
The lattice parameterization of the sequence y(k} is unique. That is. 


the lattice parameters given by (6.13) are unique. 
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b. Proof (Outline) 

The truth of this theorem is a consequence of the fact that the lattice 
can be interpreted as a Fourier series expansion of the sequence y(k) in terms of 
the orthonormal backwards errors (see equation (6.17).) The uniqueness of the 
Fourier series is well known [Ref. 47]. 

3. Synthesis Models 

A model analogous to that illustrated in Figure 4.7 can be used to 
regenerate the original sequence. This requires knowledge of the forward error 
residual sequence. The other inputs (zero order forward errors) required can all 
be obtained from past estimates of the sequence or from initial conditions. 

The probability distributions of the output forward error sequences must 
be known if a synthesis model is to be constructed which accurately reproduces 
the original signal statistics. It is not clear that any general statements can be 
made about the nature of these distributions. It has been shown |Ref. 48: p. 357] 
that in the continuous case they will be gaussian and of the same variance as the 
original additive gaussian noise source. In the same reference it is also shown 
that for the discrete case the error sequence will not be gaussian although it will 
be white. The inaccuracies introduced through the use of gaussian noise need to 
be investigated. The stability of these lattice models is also left as an open 


question. 


fees LATION RESULTS 

In order to verify the theory, a computer simulation was performed. Uniform 
white noise was used to excite the system. This choice provides a bounded input 
so that the response of the system used for the simulation could be guaranteed to 
be stable for a suitable choice of system parameters. A sufficient condition for 
the AR model to be stable. given a driving signal whose absolute value is 
bounded on (-a.a). where 0 < a < 1. is given by 

r F 


This condition guarantees that the output never exceeds unity and thus remains 
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Figure 6.1: Quadratic Nonlinear Lattice Filter 
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bounded. This condition is extremely restrictive. however, it does provide a 
simple rule for building models which can be used for simulation purposes. (The 
issue of the inaccuracy introduced by using uniform noise to drive the system has 
been avoided.) 


The model tested was 


-—.1 .22 .02 
[Hy,,] = | -02 .09 .001 (6.22) 
2 05 —.03 


Using the nonlinear form of the Levinson algorithm, equation (6.11) the 


following model parameters were estimated for two repetitions of the experiment, 


—.1003 .26038 .04568 


(Hy a, = | -01402 .10198 03877 (6.23) 
1998 .03230 —.03965 
~.1009 .2612 .06213 

[Hy a, = | 01151 .09917  .03125 (6.24) 


.19020 .04199 —.01432 


lol 


VIL CONCLUSIONS AND DISCUSSION 


A. SUMMARY OF NEW RESULTS 

In this dissertation the use of tensor concepts in the field of signal processing 
was investigated. The research was successful in a number of areas, extending 
known results and introducing some new ones. In particular, tensors were used 
to study nonlinear and multidimensional systems. 

The nonlinear modelling problem was formulated using tensors. By 
interchanging the roles played by the time order and the nonlinear order a new 
form. different from (although equivalent to) the traditional Volterra Series was 
developed. Using this new form, tensor equivalents of the discrete-time Wuiener- 
Hopf and Yule-Walker equations were derived. These equations can be solved for 
the unknown system parameters. Conditions for the solution of the Wiener-Hopf 
equations. without requiring matrix inversion. were specified. This resulted in a 
computational saving of O((p~1)*%-")) operations, where N is the largest time 
delay present and p is the highest exponent present in the system model. It was 
further shown that linear adaptive algorithms. such as LMS and RLS. can be 
applied to solve for the system parameters. 

Existing Lattice filter algorithms were reformulated in a tensor framework. 
It was shown that they can be considered to be equivalent expansions in 
alternate coordinate systems. These results were then applied to the solution of 
the 2-D autoregressive modelling problem. Several new 2-D lattice structures 
were deduced. These models are not efficient in the sense that an AR model 
possessing O(N*) parameters would require O(N‘) parameters when recast into a 
lattice form. It was shown that with proper assumptions of shift-invariance the 
complexity of the lattice inedels can be reduced to O(N*) parameters. 

The 2-D lattice models developed were then used to deduce a nonlinear 
lattice model. This model] differs from that of other researchers in that it allows 


updates to be performed in the nonlinear order as well as time order. 


Finally. it was shown that these lattice models are amenable to systolic array 


implementations. 


B. FUTURE DIRECTIONS FOR RESEARCH 

The multidimensional and nonlinear lattice theories are by no means 
complete. So far, several new lattice models have been developed. It now 
remains to investigate the properties of the models. 

For the 2-D lattice filters, conditions for model stability have vet to be 
established. The stability of multidimensional systems is a much more difficult 
problem than for the 1-D case because it is usually impossible to factor 
multivariate polynomials. However, necessary and sufficient conditions for AR 
2-D model stability have been established. It is believed that these can be used 
to derive lattice stability conditions. 

Another important topic is the synthesis of 2-D transfer functions using 
lattice structures. Stated differently, the problem is to design a 2-D lattice 
parameter filter that will implement ai given 2-D frequency response 
characteristic. 

For the nonlinear lattice structures similar issues need to be addressed. In 
this case. however. the problems are more complex since the stability and 
synthesis questions have not been solved for the nonlinear AR models. 

Before the synthesis problem for the lattice implementation can be solved 
several more basic questions need to be answered. The first is the definition of a 
useful and meaningful specification of the desired filter characteristic. The 
nonlinear AR filter affects not only the frequency content of the driving signal 
but also its probability distribution. The filter will also have different frequency 
characteristics for different driving signals. All of these effects should be included 
in the specification. 

It has been shown that for the continuous case. a nonlinear whitening filter 
will yield a signal with a gaussian probability distribution [Ref. 48: p. 357]. For 
the discrete case it has been shown that this is not true. The inaccuracy 


introduced by approximating the input to the synthesis model by a gaussian 


1538 


introduced by approximating the input to the synthesis model by a gaussian 
signal is an important question. 

Stability of the nonlinear lattices is also an important question. No simple 
stability tests exist for the nonlinear AR models. Their recursive nature makes 
stability difficult to analize. In general, stability will be a function of the input 


as well as the system which significantly complicates the problem. 
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APPENDIX A: ALTERNATE PROOF OF THEOREM 4.2 


In this appendix an alternate proof of the generalized error order update 
equations (4.40) is provided. The proof presented here relies on geometric 
arguments which are possible if an abstract mathematical framework is adopted. 
In this appendix, a random process will be considered to consist of a time-series 
of random variables. Each of these random variables is thought of as a vector in 
an infinite dimensional Hilbert space. For a rigorous discussion of these concepts 


see for example [Ref. 51,52]. 


A. DEFINITIONS AND FORMULATION 
A discrete random process Y = {rt O<i< x} can be considered to span a 


K+1 dimensional subspace S**! of an infinite dimensional Hilbert space S 
(assuming completeness and the linear independence of the y(i)). The inner 


product on this space is defined to be 


<u,v> = ef} (A.1) 


where u and v are elements of S. This inner product induces a norm 
uo} = <u.u>l/? (A.2) 


The error, e<,. in predicting the element y(k+1) from the previous N 


elements of Y 1s given by 


k 4 
ena = dy bS(k-1ytA) (4.3) 
\=0 
where 
aes 2 J : eas =) hy es hipaa ; ie [iO — "0 (A.4) 


A normalized version of the forward error is given by 


ee 
Cues 
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Ken 
= Yay) (A.db) 


where 


1 


N a 
oe 


hy (k+1) (A.6) 


The backwards prediction error, ry, is the error associated with the 


prediction of y(i) given the next N elements of Y. It is given by. 


re 
new = 3D hy (k-N)y(A) (A.7) 
a= 
where 
My (kK-N)i = [0 --- 0 1 —hyina —hy-ny2 -°° —hy 0 --- 0} (A.8) 


A normalized version of r,X,, can be defined as 


N 
Pe-N 


inauee = aaah, (A.9a) 
TKN | | 
K ’ 
= J bA(kK-N)y(2) 2) 
A=0 
where 
| h*(k—N 
Pik=N) = —— : y (A.10) 
Thon | | 


By orthogonal we mean that given two elements of the space S. u and v. say. 


then 


0 foru #v 
<u,v~ -{ (A.11) 


on wy Boru = v 
This will be indicated symbolically as 
: ae ( Asal) 


The expression 
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u_ S™ (Ams 


_——— 


implies that u is orthogonal to all elements of the subspace S™. 


The symbol ‘© is used to indicate direct sum. For example 
S=s' ©3S (A.14) 


means that the space spanned by S* is the the space spanned by the Cartesian 
product of the underlying sets of the spaces S! and S? [Ref. 47: p. 196]. 

Finally, the symbol ’v’ will be used to mean the span of. 

These definitions are sufficient to state and prove the generalized error order 


update theorem. 


B. ERROR ORDER UPDATE REGUSIO. = 
1. Theorem 4.2 
The (N+1) order errors can be calculated from the N- th order errors 


through the recursion 


k+] N+1 +1 eS) 
eral (E + | : ; = (A.15) 
Tien </ie (pau PN+1 rk_N 
where 
PN+1 = SOx 41-Ph-n? (A.16) 


is called the partial correlation coefficient or reflection factor [Ref. 17: p. 


37]. 
2a eroor 
Mamie tiat 
ae Se {si N= 1) 4 (kN 2) sth (A.17) 
but 
rox = Sy eee ples Uk eee a (.A.18) 
Therefore. 
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ot = Vv fis y(k—N- 1) ne vu (A.19) 
which can be written as 


Spt = Ss." ®v fs} (A.20) 


For the updated forward error, e,\4!, the following is true 


a | s} _» fits] (A.21) 
Since eX, , S we need only orthogonalize it (by a Gram-Schmidt procedure) 


with respect to Y ‘e ; focobpaim «,. Lhus. 


N+] 
end = C41 — CoN [A 22) 
where 


N —N 
ey Al ee N> 


s ~ (A.23a) 
rk-N 
= PN] [K-i { (A 23b) 
Therefore 
| 
# C1, Es - 
=o Ne k+1 aN k — 
et - ~~ Ney vy K-17” PN=+1 TR-N (A.24) 
CK +) 


Toon the jact that 


eye] 1 


a a a 25} 
ee a ie. 25 
hel / T= loxeal 


will not be repeated here (see Chapter 4. Section B). 
Similar arguments can be used to verify the backward error order update 


equation given in (A.15). 
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APPENDIX B: FORTRAN PROGRAM LISTINGS 


CF RE RE KR RRR RARER FH KA A AAI A A eat Nace ce 


2DLAT MAIN PROGRAM 
PROGRAM TO TEST 2D LATTICE SUBROUTINES 


WRITTEN 29 APRIL 1985 


C 
C 
C 
C 
G 
G 
CREF ERE REESE ERE EL ES ESS EH KR Ee ee ree 

REAL*8 ALPHA(26,26), BETA (26,26), A (26,26) ,B(26,26), RHO(26,26) 

REAL*8 HZ(25). ¥ (256,256),R (26,26) 

REAL*4 GAIN(40,40) 


DEFINE THE AR FILTER COEFFICIENTS 


DAT we eGre 22. 12.1112) 0.07 

DATA HZ/1.0,-.470,.007,.295,0.0,.015,.055,.003,.022,16*0.0/ 
DATA HZ 1.0,.03,-.015,-.011,.033,-.47,.195,.055,0.0,-.085,0.0.0.0, 
* 003,.022.-.0001,.0067,0.0,0.0,0.0,0.0,0.0,0.0,-.002,0.0,.0001 / 


Cy One@ 


DEFINE OTHER PARAMETERS 


MODEL SIZE 


err erenre, 


N = 
MN 


Hen 


N-N 


ACTUAL SYSTEM SIZE 


") Clg 


NA = 5 
VENA = San 


G NUMBER OF POINTS IN THE PLOT IS IP X IP. 
IP = 40 
CG NUMBER OF DATA POINTS TO ESTIMATE CORREEARIONS 150) 5). 
C 
TYS - 12& 
i CALL APPROPRIATE SU BROMIIES 
CAL). TiN 72 N ACG Ao) 
CALL PLOPTRIGALNIP TP. ORIGINAL Sree iia 1s) 
CALL GERRATGIZ NAY iy 3) 
CALL CORMIER ars iy o-RaN) 


CALL SCHUR(RHO.R.ALPHA,BET A. MIN) 
CALL LEVSON AB ROT vies) 
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10 


DO 101 = 1, MN 
HZ(I) = A(1.1)/ A(LJ) 

CONTINUE 

CALL TXFCN(HZ,N.GAIN.IP.1) 

CALL PLOT3(GAIN,IP.IP,'4-TH ORDER LATTICE APPROXIMATIONS’) 

STOP 

END 
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CFF ER EERE RARER EEE EE Be Fe ee ee 


SUBROUTINE GENRAT 

THIS SUBROUTINE GENERATES A 2D DATA FIEED FROM Tie 
FILTER COEFFICIENTS. A WHITE NOISE INPUT AND WHITE NOISE 
INITIAL CONDITIONS ARE USED. 


WRITTEN 29 APRIL 1985 


RAKRAKAKAKKAKKAAKAKKKAAKAKKAKKKKKAKAKAKKKKAKAKAKKAKAKKAAKKKKKRKKKKKKKAKKKKKAKAKAKKKAKE 


OO G1 Cie @ Geo) 


SUBROUTINE GENRAT(HZ,N,Y,IYS) 
REAL*8 HZ(25),Y (256,256), ADD 


FETCH THE RANDOM NUMBER GENERATOR SEED 
UNDER FORTHX, STORE SEED IN FILE 13. 
UNDER DISSPLA ENTER SEED EXPLICITLY. 


READ(13,10) IY 
10 FORMAT(I10) 
REWIND 13 
IY = 817346 


SE S7E7e SL eer ®, 


GENERATE THE RANDOM FIELD 


QV ene 


MINTS IN 
MA I= MIN=1 
DO 40 f=71Ys 
DOr 0r = ives 
rile) RAND (ieee 
NIE a6) 
eat 
DO 20 K = 1.MNMI]1 
IE (MOD(K IN) NEOWGOTO ms 
hie = 1 
ed 
GO TO 16 
15 he = Aa] 
16 Pte) LEO) OR (CER ee Ce © Omg 
AD POS (i-Ki.J-Kd) 
GO TO 1& 
A NOD = LRAND(TY)|=25 
No eee). 3) = TZ Ie iy 
CON PINE 
CON Tine 
CONTINGE 


CS — — 
~~ — 
— — — “YO 


STORE Tile \ DON SER See 


WRITE(13.50) TY 
0 FORMAT(I10) 


ho Ve We ae eee 
qn 
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(EEA GE RRR EER ERR OO a eee 


YO @ SO GO O1ee oe © 


Cree 


DiEOE®. 


10 


of 


SUBROUTINE GoRnuAtr 

THIS SUBROUTINE PRODUCES A CORRELATION MATRIX FROM 
A 2-D DATA FIELD IN AN ORDER 

WHICH IS COMPATIBLE WITH SUBROUTINE SCHUR 

WRITTEN 30 APRIL 1985 


fete ee eee ee ee SSeS SSS SSS SSS SS SS SET SSS SSS SESS SESS SSS SS SESE ST SS STS SS SS SS FS FS 


SUBROUTINE CORLAT (Y,IYS,R,N) 
REAL*8 Y(256,256),R(26,26),SUM 


DEFINE CONSTANTS 


MN = N*N 


BEGIN OUTER LOOP 


DO 100M ir JN 
Nid MP1l2a1 
LEIM = 2" MP1 -1 
DO 90 L = 1,LLIM 


Fo = - 1 

I1 = MO 

Jie 9072 

1 (MOD/(LO,2).EQ.0) GO TO 10 

| ae | 

EN 10 

JRe= IR. 

IR=IR-+ 1 

KEBOT = L 

DOesGe el = \IPTN 
Wea NP] = t 


KUIMisge2” NPI - 1 
DbO 70K = KBOT.KLIM 


iO: = K-! 
i 
j= K0; 2 


IF (MOD(K0.2}.EQ.0) GO TO 20 
ee. 


eae NG 
Leo — 11 
ONS eel 


KIMAX =TS 

Kivi = 1OFF = 1 
iblOkr GT.0) GO TO 20 
KIMAX = [YS + IOFF 
Kot iN =" 1 

K2nrAy = ys 
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K2MIN = JOFF + 1 
IF (JOFF.GT.0) GO TO 40 
K2MAX = IYS — JOFF 


K2MIN = 1 
40 SM = 0:0 

JR = JR-+1 
G WRITE(6,44)IR,JR,I1,J1,1,J, JOFF,JOFF 
C44 FORMAT (8(2X,14)) 


DO 60 K1 = K1MIN,KIMAX 
DO 50 K2 = K2MIN,K2MAX | 
SUM = SUM + Y(K1,K2)*Y(K1-IOFF,K2-JOFF) 
50 CONTINUE 
60 CONTINUE 
R(IR,JR) = SUM/(KIMAX-K1MIN+1) /(K2MA X-K2MIN-~ 1) 
70 CONTINUE 
KBOT = 1 

80 CONTINUE 
90 CONTINUE 
100 CONTINUE 
C 
C FILL IN THE SYMMETRIC HALF OF CORRELATION MATRIX 
G 

DO 1201 = 2.MN 


IM1=I1-1 
DO 110 J = 1,IM1 
R(I.J) = R(J,I) 


110 CONTINUE 
120 CONTINUE 
c 
G 
RETURN 
END 
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Cl EER RRR ERR RAAAAE RRR RE ee eee ee ee ee 


SUBROUTINE TXFCN 


GENERATES A 2-D FREQUENCY RESPONSE GIVEN THE TRANSFER 
FUNCTION COEFFICIENTS 


WRITTEN BY DR. B. MADAN 
MODIFIED BY P.J. LENK 29 APRIL 1985 


PARAMETER KO INDICATES THE ORDERING OF THE COEFFICIENTS 
- KO = 0: PARAMETERS IN ROW-MAJOR ORDER 
- KO = 1: PARAMETERS IN TWIDDLED ORDER 


ee oe eo oo oo oo eK RK KK KK KK KK 


OCY OO O1@.O GeO -G oe Gre 


SUBROUTINE TXFCN(HZ,N,GAIN.IP,KO) 
REAL*8 HZ(25),CONVRT(5,5) 

REAL*4 GAIN(40,40) 

COMPLEX CSUM 

COMPLEX*8 ARG1,ARG2 


C. DEFPINEZCONSTANTS 
PI = 3.14159 


DETERMINE ORDER AND REORDER IF NECESSARY THE COEFFICIENTS 


ere) 


IF (KO.EQ.0) GO TO 60 
pee — 1), 
DO°26 ME I= 1X 
MO = MP1-1 
LES" MP i= 1 
DOs2Z0t= Luin 
LOR JL = 1 
LP eNi0 
bebo LG a 
IF (MOD(LO0.2).EO 0) GO TO 16 
Le 
J=a10 
10 JR Sk J 
CONN 1.J—1) — HZeie) 
Dt) CON TIE 
a0 6 oC CON PI 
C 
G TRANSFERS CORFPICIENT ob ee loa 
% 
Jit aee 
DO 50 1 25a 
DO 407 ie 
JR =n 
HZ( JR eONy RT (Ty) 
WRITE Ge. 123) 


35 FORMAT(2X.E12.5) 
40 CONTINUE 
50 CONTINUE 


C 
© PROCEED WITH TRANSFER FUNCTION EVALUATION 
C 


60 DW = 2.0*PI/FLOAT(IP-1) 
DO 100 IW1 = 1,IP 
W1 = DW * (IW1 - 1)- PI 
DO 90 JW1 = 1,IP 
Al = Wil 
W2 = DW * (JW1 - 1)- PI 
17 = 0 
CSUM = CMPLX(0.0,0.0) 
Os80.1— TN 
ARG1 = CMPLX(0.0,-A1?) 
Al = Al + W1 
A2 = W2 
DO 70 J = 1,N 
IZ=1Z+1 
ARG2 = CMPLX(0.0,-A2) 
A2 = A2 + W2 
cS WRITE(6,77)1,J,ARG1,ARG2,IZ,HZ(1IZ),CSUM 
CTi FORMAT (2(2X,13),4(2X,E£12.5),2X,13,3(2X,E12.5)) 
CSUM = CSUM + CMPLX(SNGL(HZ(IZ)),0.0)*CEXP(ARG1+ARG2) 
70 ~CONTINEE 
80 CONTINUE 
GAIN(IW1.JW1) = 1.0;CABS(CSUM) 
GAIN(IW1.JW1) = GAIN(IW1,JW1)*GAIN(IW1,JW1) 
c WRITE(6.78)IW1.JW1.CSUM.G AIN(IW1,JW1) 
C78 mete veal (2t25-13).3(2\.E12.5)) 
90 GONTINUE 
C mr Ve, 1e)IW 1 (GAIN(TW 1,1) J=1,IP) 
C18 BORAT x 13 75| 2X .E12.5)) 
100 CONTINUE 
Fest tURN 
PND 
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SUBROUTINE PLOTS 


THIS IS A ROUTINE TO USE THE DISSPLA PACKAGE TO 
DRAW A THREE DIMENSIONAL PLOT OF A 2-D FILTER’S 
FREQUENCY RESPONSE. 


WRITTEN BY DR. B. MADAN 
MODIFIED BY P.J. LENK 29 APRIL 1985 


SPCC SSeS SS SSeS eS SSeS SSS SSS SES TESS SESS ESSE ESS SSS SES SSS SSS SEES SSS SES SF 
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SUBROUTINE PLOT3(A,IM,JN,LABEL) 
DIMENSION A(IM,JN) 

INTEGER LABEL(100) 

CALL TEK618 


INITIACIZE ee PLOTTING DEVICE 


CaiGee) 


WRITE(6,51) 
FORMAT(! 1 CALL TEK618 OK’) 
CALL THE DEVICE 
CALL RESET(’ALL’) 
WRITE(6,52) 
FORMAT(’ 2 RESET ALL OK’) 
CALL SETUP FOR CUBE 
Al=FLOAT(IM) 
A2=FLOAT(JN) 
CALL PAGE(11..9.5) 
W RITE(6,53 
53 FORMAT(’ 3 CALL PAGE OK’) 
CG see Ghat >) 
C CALL PHYSOR(0.,0.) 
CALL AREA2D(7.0.7.0) 
WRITE(6.54) 
54 FORMAT( 4 CALL AREA 2-D OK’) 
CALL SIMPLX 
CALL HEIGHT(.2) 
CALL HEADIN(LABEL.100.4.1) 
CA eietenaea. 14) 
CALL VIEW (-10.0.-5.0.15.0) 
CALL VOLM83D(12.0.12.0.12.0) 


si 


ae 


( CXS VAI LABELLING RO lias 
CMB ee Nr 

C CVE eV NIS LABELING 17 iE 
eh. la aaNet: =a 3 

G Clezes IS LABELLING KOU Ri 


CALL Ayal 22) 
Chere SCRE ACE VLOT Ket tie 


= 


CALL GRA 3B mun 20. 1.0.0. 2.180.010 21aulen) 
1500 -CAELE Sie ee reel 0) 
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SUBROUTINE SCHUR 
CALCULATES THE REFLECTION FACTORS FROM THE CORRELATION MATRIX 


WRITTEN 29 APRIL 1985 
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SUBROUTINE SCHUR(RHO,R,ALPHA,BETA,N) 
REAL*8 RHO(26,26),R(26,26), ALPHA (26,26), BET A(26,26), RNORM,T 


Cc 
C INITIALIZE THE ALPHA AND BETA ARRAYS 
cc . 
DO 101=1,N 
DO5J=1,N 
ALPHA(I.J) = R(I,J)/DSQRT(R(I,])) 
BETA(l.J) = ALPHA(I.J) 
RHO(I.J) = 0.0 
5 CONMINUE 
é W RITE(16,7)(ALPHA(I,J),J=1,N) 
Gn FORMAT(5(2X.E12.5)) 
10 CONIINGE 
C 
C BEGIN CALCULATING THE REFLECTION FACTORS 
C 
DO 50 J = 2.N 
NJL=N-J+1 
DO 401 = 1.NJ1 
Mi 4124 
IPl=1-1 
RHO(I.JI1) = ALPHA(I,JI1)/BETA(IP1.J11) 
RNORM = DSQRT(1.0 - RHO(I,J!1)*RHO(I,JI1)) 
DO 30K = 1.N 
T = ALPHA(LK) 
ALPHA(I.K) = (ALPHA(I.K)-RHO(L.JI1)*BETA(IP1.K)), RNORM 
BETA(I.K) = (BETA(IP1.K)-RHO(I.JI1)*T) ‘RNORM 
30 CONTINUE 
10 © CON 
C WRITE G.42)5.((ALPHA(LK).K—1.N).E- 1) 
C42 ORIN 2X.19.4(2N.E12.5)) 
C WRITE(16.42)J.((BET A(1.K).K-=1.N).1=1.N) 


AQ CONSTI Gale 
( 


RET ines 
Je Sle: 
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SUBROUTINE LEVINSON 


GENERATES THE AUTOREGRESSIVE MODEL COEFFICIENTS GIVEN THE 
REFLECTION FACTORS 


WRITTEN 29 APRIL 1985 


-f eee SEES SESE SESS SESS SSS ESSE SSS SSS TES SSE SESE SSE SESE SESE SSE SSS STS SEE SESS ST SS SE tS 
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SUBROUTINE LEVSON (A,B,RHO,R,N) 
REAL*8 A(26,26),B(26,26),RHO(26,26),R(26,26),RNORM,T 


C 
C INITIALIZE THE A AND B MATRICES 


C 
DO 201 = 1,N 
DO 10J =1,N 
A(1,J) = 0.0 
B(I,J) = 0.0 


10 CONTINUE 
A(LI) = 1.0/DSQRT(R(L])) 
B(I,]) = A(I,]) 

20 CONTINUE 


C 
C CALCULATE THE AR PARAMETERS USING LEVINSON’S ALGORITHM 
C 


wo 60 J = 2,N 
iNol= N-jJ+4+ 1 
wo 50] = NII 
fie 1=— jJ- ] 
RNORM = DSQRT(1.0 - RHO(I,IJ1)*RHO(I,1J1)) 
WO 401K =1.1351 
T= A(loky 
Perm — (ALK RHO(lIS1)*B(I41.4))/RNORM 
Bie = (B(l- 1K) - RHO(I.1J1) SB) “RNORM 


40 Contin UE 

mie 4=©6CONTINUE 

C WRITE(16.77)J.((A(IL.K).K=1.N).J=1,N) 
c RT aiees 0)d.((Bil,.K).K=1.N).1=188) 


Cis FORMAT(/2X.12.4(2X.E12.5)) 
60 CONTINUE 


mo 66] = 1.N 
eR NU = Afi.l) A().1) 

(- Pvt. oo | R NORA 
C WRITE(17.65)RNORM 
C65 PORMAT(24.§E12.5) 
a ONTINUCE 
(: 
C 

RETURN 

END 
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FUNCTION URAND 
TAKEN FROM "COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS" BY 


G.E. FORSYTHE. M.A. MALCOLM, AND C.B. MOLER 
PRENTICE-HALL, ENGLEWOOD CLIFFS, NJ., 1977, P. 246. 


KKKKKKKKKKKKKKKKKKK KKK KKK KKK KKK KKK AK KKK KKK KKKK KK KKK KKK KKKKKK KKK KA KKK KK 
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REAL FUNCTION URAND (LY) 
INTEGER IA,IC,ITWO,M2,M,MIC 


C  URAND IS A UNIFORM RANDOM NUMBER GENERATOR BASED ON THEORY AND 
C SUGGESTIONS GIVEN IN D.E. KNUTH (1969), VOL. 2. THE INTEGER IY 
C SHOULD BE INITIALIZED TO AN ARBITRARY INTEGER PRIOR TO THE FIRST 
C CALL TO URAND. THE CALLING PROGRAM SHOULD NOT ALTER THE VALUE OF IY 
C BETWEEN SUBSEQUENT CALLS TO URAND. VALUES OF URAND WILL BE RETURNED 
C IN THE INTERVAL (0,1). 
C 

DOUBLE PRECISION HALFM 

REAL S 

DOUBLE PRECISION DATAN,DSQRT 

DATA M2/0/,ITWO/2/ 

IF(M2.NE.0) GO TO 20 
Cc 
C IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORD LENGTH 
C 

M=1 

10 M2=M 

M=ITWO*M2 

IF(M.GT.M2) GO TO 10 

HALFM=M2 
Ec 
C COMPUTE MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD 
G 


[A=8*IDINT(HALFM*DATAN(1.D0) /8.D0)—5 
IC =2*IDINT(HALFM?(.5D0-DSQRT(3.D0)/6.D0))+1 
MIC =(M2-IC})— M2 
C SIS MME SCALE FACTOR FOR CONVERTING 1 oO Flos TIN CcHren a 
S=.5 HALFM 


G COMPU ENT BRANDON Bik 


20 1) Ses 
z 
C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHICH DO NOT ALLOW 
C INTEGER OVERFLOW ON ADDITION 
G 
IF(1Y.GT.MIC) 1Y=(IY-M2)-M2 
C 


Iy=1¥+ICc 


CAO ea 


THE FOLLOWING IS FOR COMPUTERS FOR WHICH THE WORD LENGTH 
FOR ADDITION IS GREATER THAN FOR MULTIPLICATION 


IF (TY /2.GT.M2)ITY=(TY-M2)}-M2 


THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER OVERFLOW 
AFFECTS THE SIGN BIT 


IF(IY.LT.0)TY=(1Y +M2)+M2 
URAND=FLOAT(IY)*S 
RETURN 

END 
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PROGRAM NLMAIN 
THIS IS THE PROGRAM TO TEST THE NONLINEAR LATTICE MODEL 


WRITTEN 7 MAY 1985 


KK KKK KKK KK KK KKK KKK KKK KK KK KKK KKK KKK KKK KKK RK RK KKK RK KK KK KK KKK KKK KKK KKK 


REAL*8 ALPHA(26,26), BET A (26,26), A (26,26), B(26,26),RHO(26,26) 
REAL*8 HZ(5,5), Y(10000),R (26,26) 


DEFINE SYSTEM PARAMETERS 


DATA HZ/.2..02,-.1.0.0,0.0,.05,0.09,.22,0.0,0.0,-.03,.001,0.02,0.0 


* 0.0,10*0.0/ 


DEFINE SYSTEM CONSTANTS 


XO A ae 
K =6-] 
WRITE(6,4)(HZ(K,J),J=1,5) 
FORMAT(5(2X,E12.5)) 


CONTINUE 
T¥S = 1000 
NA = 3 


MINA = NA * NA 
MNAP1 = MNA + 1 


DEFINE MODEL CONSTANTS 


a 
MN=N?N 
MNP] = MN + 1 


CALIFSUBROUTINES 


CALL NLGEN(Y.N.HZ.IYS) 
CATION EOI dh V1YS Ra) 
DO 301 1.MNPI 
W RITE(16.20)(R(1.J).J=1,.MNP1) 
FORMAT(:.5(2X.E12.5)} 
CONMANT E 
CALL SCHUR(RHO.R.ALPHA.BETA.MNP1) 
DO 60 I = 1.MNP1 
W RITE(16.20)(RHO(I.J),J=1.MNP1) 
CONTINUE 
CALL LEVSON(A,B.RHO,R,MNP1) 
DO 70] = 1,MNP1 
W RITE(16,20)(A(I.J),J=1,MNP1) 
CONTINIS 
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SUBROUTINE NEGEN 


THIS SUBROUTINE GENERATES AN OUTPUT SEQUENCE FROM THE SYSTEM 
DESCRIBED BY THE MODEL PARAMETERS CONTAINED IN H(,). IT 

USES WHITE NOISE UNIFORM ON (-.5,.5) TO EXCITE THE SYSTEM. 

THE INITIAL CONDITIONS ARE ALSO DRAWN FROM THIS DISTRIBUTION. 


WRITTEN MAY 7 1985 


(tet ec eee S SPSS SS SSS SSS SSS TESS SSS SSS SSS SSS SSS SSE SS TESS SS SSS SS SS TS 


OO: @i@  O1@ 2:6, a 


SUBROUTINE NLGEN(Y,N,H,IYS) 
REAL*8 Y(10000),H(5,5),DSEED 


C 
C FETCH THE RANDOM SEED 
c 

READ(13,10) IY 


10 FORMAT{(II0) 
REWIND 13 


Cc DSEEDE— WD PEOATL (TY) 
C 
C SET UP THE INITIAL CONDITIONS 
C 

Y (fee CURAND 5) 

Y (22. (URAND (DY) 5) 
C CALL GGNML(DSEED,JYS,Y) 
S: 
C CALCULATE THE REMAINING VALUES OF THE SEQUENCE 
© 

DOaOl= silys 
CG Y= 2. (RAND es) 
C ye — ct) 

Yo — URANDITY) 

Cc DOraC sla i \ 
(© US ek | 
(o DOO kh Ie 
ce Kove -1 
@ YI] Yl) - H{i2.K)*COORD( Yee Isl) COORDS (1-2) ham, 
C20) CONTIN 
C30 CONTA E 


Q CONTIN 


PL Nissi 


Oe ee Ne 
rr er | 


(y= DIN Seb) 
WRITE(13.50)1Y 

50 FORMAT{IIO) 
REWIND 13 
RETURN 

END 
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SUBROUTINE NLCLAT 

THIS SUBROUTINE PRODUCES A CORRELATION MATRIX FROM NONLINEAR 
TIME SEQUENCE IN AN ORDER WHICH IS COMPATIBLE WITH SUBROUTINE 
SCHUR. 


WRITTEN 7 MAY 1985 
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SUBROUTINE NLCLAT (Y,IYS,R,N} 
REAL*8 Y(10000),R(26,26),SUM,VEC(26) 


DEFINE CONSTANTS 


wero oO) 


MN = N*N 

MNP1 = MN + 1 

Iyeeni2 = LYS - 2 

FIYSM2 = FLOAT(TYSM2) 


INITIALIZE R MATRIX TO ZERO 


oO @ 


DO 201 = 1.MNP1 
DO 10 J = 1.MNP1 
R(I,J) = 0.0 
10 CONTINUE 
20 CONTINUE 


ee) BEGIN OUTER LOOP 


merso | = 3,1YS 
i 1 
mEC(IR) = Y(I) 
WO sOe1P1] = 1.N 
MO = MP1- 1 
ire 29NiPi - 1 
PO 4 = 1 L LIM 
Poe =| 
eee ait) 
Jie" i0)).2 
PeronotEt 2) EO.0) GO TO 30 
|) Vere 
i AIC 
30 IR = IR - 1 
VEC(IR) = COORD(Y(I-1).11)*COOR D(Y(I-2}.J1) 
40 COIN E 
) CONTINUE 


o( 
C 
© CALCULATE THE CORRELATIONS 
C 


Om) = TAN 
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DO 60 K = J,MNP1 
R(J.K) = R(J,K) + VEC(J)*VEC(K) 


C WRITE(6.12)VEC(J), VEC(K).R(J.K) 
OMe FORMAT(3(2X,E12.5)) 
60 CONTINUE 


70 CONTINUE 
80 CONTINUE 
C 
C DIVIDE BY THE NUMBER OF DATA ELEMENTS CONSIDERED 
Cc 

DO 100 J = 1,MNP1 

DO 90 K = J,MNP1 
R(J,K) = R(J,K)/FIYSM2 

90 CONTINUE 
100 CONTINUE 
C 
C FILL IN THE SYMMETRIC HALF OF CORRELATION MATRIX 
€ 

DO 1201 = 2,MNP1 


Mie ed 
DO MGH .aieive 
R(L.J) = R(T) 


110 CONTINUE 
120° COMIN UE 
C 
C 
RETURN 
END 
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NONLINEAR WIENER MODELLING PROGRAM 
iiieeeor> FUNCTIONS CONTAINED IY ROUTINE COORD 
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C 
C 
C 
C MODIFIED 1 JAN 85, 14 JAN 85 
C 
C 


DOUBLE PRECISION A(25,26),Z(25),X(15000),Y(15000), VAR,XA,YA,YHAT 
DOUBLE PRECISION ZZ(25),XM1,SEED,STORE(25) 


INPUT SIMULATION PARAMETERS 


yore 


READ(13,77)TY 
REWIND 13 
WRITE(6,40) 
40 FORMAT(2X,’MAGNITUDE OF NOISE’) 
READ(5,41)VAR 
41 FORMAT(F12.5) 
WRITE(6,42) 
42 FORMAT(2X,NUMBER OF POINTS’) 
READ(5,43)N 
43 FORMAT(I5) 
WRITE(6,44) 
44. FORMAT(2X,,MAXIMUM POWER OF X ’) 
READ(5,45)LW 
45 FORMAT(I1) 
WRITE(12,46)N 
WRITE(6.46)N 
for oR MA T(2A THE NUMBER OF POINTS WAS ’,I5) 
WRITE(12.47)LW 
WRITE(6.47)LW 
47 FORMAT(2X’°THE MAXIMUM POWER OF X JIS ’,I1) 
WRITE(6.48)VAR.VAR 
PORITE(12.48)VAR.VAR 
fae ORMAT(2X” THE NOISE IS UNIFORM FROM -’,F9.5,’ TO +’,F9.5) 
Bay = LW — 1] 
Peak = VAR“*2.0 
Mois 2 | = 1.25 
Sl ORE(I) = 0.0 
fez CONTINUE 
mai} = VAR*(URAND(IY) - .5) 
med) = 1.0 
WO 2 1=2.N 
il) = Vak (Ul RAND(IY) - .5) 
my = UNKNO WM (X(I).N(I-1}} 
Se OO NTINUE 
Pale NCRLAT(A.Y.A.LSOW,LW.N) 
C WRITE(6,134) 
C WRITE(12,134) 
ea FORMAT(/2X”NO INVERSION SOLUTION’) 
moO 1421-1. LSQOW 
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LL 3) = Al ESON ae) 
STORE(J) = STORE( J 273) / 25-6 
142 CONTINUE 
C DOv-Ise J = 
CG IMIN = (J-1)*LW + 1 
Cc IMAX = J*LW 
G WRITE(6,173) (ZZ(I),] = IMIN, MAX) 
c WRITE(12,173) (ZZ(I),] = IMIN, IMA X) 
C173 FORMAT(5(2X,E12.5)) 
C133) CONTINUE 
CG CALL SOLVE(A,Z,LSQW) 
\G WRITE(6,135) 
C WRITE(12,135} 
C135 FORMAT(/2X,’,USING FULL MATRIX INVERSION’) 
C DO 2732 iw 
c IMNe= ee Ww - 41 
C IMAX = J*LW 
c WRITE(6,11) (Z(I),I = IMIN,IMAX) 
Cc WRITE(12,11) (Z(1),1 = IMIN,IMAX) 
1] FORMAT (5(2X,E12.5)) 
C27 CONTINUE 
131 CONTINUE 
WRITE(6.201) 
WRITE(12,201) 
201 FORMAT(/2X"NO INVERSION SOLUTION AVERAGED OVER 25 RUNS’) 
DO 200 J-= 1.EW 
IMIN = (J -1)*LW + 1] 
[NEA = 
WRITE(6.11) (STORE(I).I = IMIN,IMAX) 
WRITE(12.11) (STORE(]).1 = IMIN-IMAX) 
200 CONTINUE 
WRITE(6,22) 
VW Roar 2 22) 
Nel "VAR (EC RANDY eo) 
DO ie = a0 
NA Bee SO 
Y A= EN KNOW [A ANM)) 
Youd 6 
VYHAT =.0.0 
DOM 1a 
BS, iene 
Nee OOn DIX MIJN) 
DO 157k 1.LW 
KA | 
No = COLD IRI) 
( WN G52) N INS 
C5, FO wes | (2(2 5. 2a i 
YHAT = YVHWAT =~ NX> N28 eZ (Jee 
YYHAT = YYHeaAT = AI” X27 ZZ ( (031) ats) 
15 CON TINIEE 
14 CON Ties 
ERROR = (YA - YHAT)‘YA 
LERKRORS 45 A=) yA Ga 
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XM1 = XA 
WRITE(6,17) YA,YHAT, YYHAT,ERROR.ZERROR 
FORMAT(5(2X,E12.5)) 
FORMAT(/6X,’YA’,12X,’YHAT’10X,, YYHAT’,9X."ERROR’,9X,’ZERROR’) 
WRITE(12,17) YA.YHAT, YYHAT.ERROR,ZERROR 

CONTINUE 

WRITE(13,77)1Y 

FORMAT(110) 

STOP 

END 
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SUBROUTINE UNKNOW 


THIS ROUTINE DEFINES THE UNKNOWN SYSTEM 


2 RK oR RK oR RRR RR ROO OR OR OE RO OR OO ROR OO RK RO RR OR RE eR ROR OR oR OR OR ROR ROK OR OOK OK KO RO OK KK OK 


COO]. G ae 


DOUBLE PRECISION FUNCTION UNKNOW(X,X1) 
DOUBLE PRECISION H(25), YHAT,XT1,XT2 


LWe=5 
H(t a2 
H(2) = -.4 
H(3)7= 208 
H(4)) = -.7 
Hi 5) = 1020 
H(Gir= 5 
H(7)= 435 
H(8) =i 
H(9) = .9 
Hi 1G) sar0 
Hii} sou 
H(12) = 1.3 
Hits an 
Hit ae 
Hilo) 20.0 
H(16) = .43 
H(i) = 
H(18) = -.05 
H(19) = .4 
H(20) = 0.0 
Si OR 
Hiezie— a. 
Hizs= 020 
Aa 20 
HiZo)e— 0 
VYHAT=00 
DOA ra Ay 
ATSB een) a 


XT2 = COORD{X1.JM1) 
DO 45K = 1,LW 
KM1l=K-1 
Mie eC ORDIN hI 
( MRM G52 1112 
RE FORMAT (2(2X.E12 5)) 
YHADMEYHAT - NT1° \T2° GRAAL = 4h} 
lh CONTINUE 
14 CONTINGE 
UNK NO emeiar 
RETURN 
END 
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FUNCTION COORD 
GENERATES OUTPUT OF THE FUNCTIONS BEING USED AS COORDINATES 


CREATED 23 AUG 84 
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DOUBLE PRECISION FUNCTION COORD (X,]) 
USE LEGENDRE POLYNOMIALS 


IF (I.NE.0) GO TO 1 
Y = 1.0 
COMO 30 
IF (I.NE.1) GO TO 2 
Y = 1.732051*X 
GO TO 30 
IF (I.NE.2) GO TO 3 
Y = 3.354102*(X*X - 1./3.) 
GO TO 30 
Y = 6.61438*(X*X*X - 3./5.*X) 


—_ 


La) 


ceo 


USE SIMPLE POWER SERIES TYPE POLINOMIALS 


Cy CPO Ory Oar Oe) a rere ere oe ar ane 


Y = 1.0 
IF (I.EQ.0) GO TO 30 
Y = X**I 

30 COORD = DBLE(Y) 
RETURN 
END 
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SUBROUTINE SOLVE 


SUBROUTINE TO SOLVE A SYSTEM OF LINEAR EQUATIONS 
USES GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING 


sk 2K KE OK ok KKK ok KK OK KK KK KK ok ok ok kK KKK OK KKK KKK KKK KKK KKK KKK KKK KK KKK KKK 


PR OROT BUG u ee 2 


SUBROUTINE SOLVE(A,Z,K) 
DOUBLE PRECISION A(25,26),Z(25),Y, TEMP,LARGE 


KM1 = K- 1 
Kei Kt 


DO 10 L = 1 KAM 
LPi =’ bal 
DOW ="LTPIK 
DARGE = DABS(A(L,D)) 
LROW =—-b 
POn2 vie LPI kK 
IF (DABS(A(M,L)).LE.LARGE) GO TO 12 
LARGE = DABS(A(M,L)) 
LROW =M 
12 CONTINUE 
C 


IF (L.EQ.LROW) GO TO 13 
DO 14M = L,KP1 
TEMP = A(L.M) 
A(L.M) = A(LROW.M} 
A(LROW.M) = TEMP 
14 CONTINUE 


2 Y = A(I.L)/A(L.L) 


DO 15 M = 1.KPI 
A(I.M) = A(I,M) - A(L.M)*¥ 
15 CONTINUE 
11 COMMUNIGE 
10 CONTINUE 


Z(K) = A(K.KP1). A(K.K} 
DO tent 
eer 
Z(P) — A(L.KP1) 
K Mie K - I 
DO17M.o1.KMI 
C= eae 
Z(I) = Z(1) - A(i,J)*Z{J} 
17 CONTINUE 
Z{1) = Z(1) A(L.D) 
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NONLINEAR CORRELATION MATRIX CALCULATOR 


CREATED 23 AUG 84 
USES FUNCTION UNKNOW TO DETERMINE FUNCTIONS FOR EXPANSION. 
X - INPUT VECTOR 
Y - OUTPUT VECTOR 
PHI - CORRELATION MATRIX 
- LAST COLUMN CONTAINS INPUT/OUTPUT CROSS-CORRELATION 


2K oe oR oR OK OK OR OK oe oe oR OR oe ok OK oe oe OR OK OE OR ok oR OE OE OR OR oR OR oe OK ok OK OK oo OKO OK OOK OK OK OK oe OK KR ROR OK KOK OK KKK KKK 


) OO GL@ @ OG Geers 


SUBROUTINE NCRLAT(X,Y,PHI,LSQW,LW,N) 
DOUBLE PRECISION PHI(25,26),X(1).Y(1),TX(5,5) 


NML=N-1 
LSQW =LW* LW 
LSQWP1 = LSQW +1 


DO 40 1=1.LSQW 
DO 39 J=1,LSQWPI 
PHI(I,J) = 0.0 
39 CONTINUE 
40 CONTINUE 


© 
6 
DO 2 i=2 SN 
NIM = Niet 
DO 3 Evi 
IM? == 1 
X2 = COORD(X(MM1),IM1} 
DO Ja rL 
JMi=J-1 
Nie COORD (XIAN ty) 
& Vera Ors ) N 1A 
C'S] POR MAMI 2(2%. 21255) 
Cee NI 2 
C Wore E680) S(K) ACIS LIT ed) 
C&O FORMAT (2X (2(2N E125), 2(2 Na Blo) 
4 COVEINnE 
2 CON Ter 
(: 
C 


BOA le 11k 
bet =i 
K See = Il 
PHI(K.LSQWP1) = PHI{K.LSQWP]) — ¥(M)*TX(LUW 
DO 7 J=1.LW 
DO 8 JJ=1,LW 
KK = (3-1)? DW = JJ 
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Bil keonen) = PHI(K KK) + TX(1,01)*TX(J,J5) 


8 CONTINUE 
7 CONTINUE 

6 CONTINUE 

5 CONTINUE 

2 CONTINUE 

C 


DO 31 I=1,LSQW 
DO 32 J=1,LSQWP1 
PHI(I,J) = PHI(I,J)/FLOAT(NM1) 
32 CONTINUE 
31 CONTINUE 


é. 

C 

C DO 171=1,LSQW 

€ WRITE(11,16) (PHI(I,J), J = 1,LSQW) 


C16  FORMAT(/,9(2X,E12.5)) 
C17 CONTINUE 
WRITE(11,19) (PHI(I,LSQWP1),] = 1.LSQW) 
19 FORMAT(//,9(2X,E12.5)) 
C 
G 
RETURN 
END 
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C ¥ x ER REA ERE RH EE a ee ete Roel ee 


NONLINEAR WIENER MODELLING PROGRAM 
THIS USES THE LMS ADAPTIVE ALGORITHM 14 FEB 84 (VALENTINE’S) 


UNKNOWN SYSTEM DEFINED IN FUNCTION UNKNOW 


KK KK AK KKK KKK KKK AK KK KKK KK KKK KK KKK KK KKK KKK KK KK KKK KKK KAKA KKK KKK KK KK KK KKK 


OOO GO] Cae 


DOUBLE PRECISION TX(5,5),H(5,5),VAR,X,XM1,Y, YHAT,SEED 
READ(13,77)1Y 
REWIND 13 
WRITE(6,40) 

40 FORMAT(2X,’MAGNITUDE OF NOISE’) 
READ(5,41) VAR 

41 FORMAT(F12.5) 
WRITE(6,42) 

49 FORMAT(2X,,NUMBER OF POINTS’) 
READ(5,43)N 

43. FORMAT(I5) 
WRITE(6,44) 

44. FORMAT(2X..MAXIMUM POWER OF X ’) 
READ(5,45)LW 

45 FORMAT(I1) 
WRITE(6.7) 

7 FORMAT(2X,,;CONVERGENCE FACTOR’) 
READ(6.8)U 

8 FORMAT(F12.5} 
WRITE(12,46)N 
WRITE(6.46)N 

46 FORMAT(2X”°THE NUMBER OF POINTS WAS °.15) 
WRITE(1247 LW 
WRITE(6.47)LW 

47 FORMAT(2X,’THE MAXIMUM POWER OF X IS °\11) 
WRITE(6.48) VAR.VAR 
WRITE(12.48)VAR.VAR 

48 FORMAT(2X.°THE NOISE IS UNIFORM FROM -',F9.5." TO +’,F9.5) 
WRITE(6.49)U 
WRITE(12.49)U 

149 FORMAT(2N “THE CONVERGENCE PACTOR WAS ’.E12.5) 
LY) =e 
VAR = VAR"2.0 
Nhl = Wong nN Diy jee oi 


(INIT EA RIZE Se PEN SOr 
DOr Ga Hoty 
DOYJ 1LW 
Hy Lee0) 
9 CON TINGE 
10. CONTINWE 
CG 


© BEGIN THETTERATION LOOr 
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C 
DO 2 K=2,N 
PhS SV AR (URAND(IY) - .5) 
Y = UNKNOW(X.XM1)} 


C WRITE(6.97).X.Y 

C97. FORMAT(2X.15.2X,E12.5,2X.E12.5) 
YHAT = 0.0 

8 


C CALCULATE THE MEASUREMENT TENSOR AND THE OUPUT ESTIMATE 
C 
DO 141 = 1,LW 
IMi =1-1 
X2 = COORD(XM1,IM1) 
DO 15 J =1,LW 


JMi=J-1 

X1 = COORD(X,JM1) 
Cc WRITE(6,52)X1,X2 
C52 FORMAT(2(2X,E12.5)) 


Tied — Soleo 
YHAT = YHAT + TX(I,J) * H(1,J) 


15 COMEENUE 

14 CONTINUE 

C 

Se CALCUPATE THE NEWVALUE OF THE TENSOR 
C 


ERROR = Y- YHAT 
DO51=1.LW 
DOm J — ile 
H(I,J) = H(1,J) + 2.0*U*ERROR*TX(I,J) 


4 CONTINUE 
5 CONTINUE 
XMi= xX 
C 
C PRINT THE NEW VALUE OF THE TENSOR 
oy 


KM1=K-1 
IF (KM1.NE.(KM1 100)*100)GO TO 2 
WRITE(6.13)KM] 
WRITE(12.13)KM1 
12 FORMAT( 2X-ITERATION NUMBER °.15) 
WRITE(6 135) 
WRITE(12.135) 
foe FORMAT(2N. THE \EW VALUE OF THE H TENSOR IS’) 
ic a 
WRITE(6.11) (H(L.J).J = 1.1) 
WRG bil2.?1) (H{1.J}J - 1.1.) 
1 FORMAT(5(2X.E12.5)) 
ay CONTINUE 
WRITE(6.22) 
WRITE(12.22} 
WRITE(6.17) Y.YHAT.ERROR 
17. FORMAT(5(2X.E12.5)) 
22 FORMAT(6X.’Y'.12X YHAT‘10X. ERROR’) 
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i) 


—s 
~1] 


WRITE(12,17) Y, YHAT,ERROR 
CONTINUE 
WRITE(13,77)1Y 
FORMAT(I10) 
STOP 
END 
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